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1  Introduction 

Screen/film  mammography  is  widely  recognized  as  being  the  only  effective  imaging 
modality  for  the  early  detection  of  breast  cancer  in  asymptomatic  women  [1].  Screening 
asymptomatic  women  using  screen/film  mammography  has  been  shown  to  significantly 
reduce  breast  cancer  mortality  [2]. 

Breast  cancer  currently  accounts  for  32  %  of  cancer  incidence  and  18  %  of  cancer 
mortality  for  women  in  the  United  States.  There  were  182,000  new  cases  of  breast  cancer 
in  the  United  States  in  1993  and  46,000  deaths.  Five  year  survival  rates  are  generally  very 
high  (93  %)  for  breast  cancer  staged  as  being  localized,  falling  to  72  %  for  regional  disease 
and  only  18  %  for  distant  disease  [3].  The  early  detection  of  breast  cancer  is  clearly  a  key 
ingredient  of  any  strategy  designed  to  reduce  breast  cancer  mortality. 

Major  advances  in  screen/film  mammography  have  occurred  over  the  past  decade  [4] 
which  have  resulted  in  significant  improvements  in  image  resolution  and  film  contrast.  Of 
major  importance  is  that  these  improvements  have  been  achieved  at  reduced  radiation 
doses.  Despite  these  advances,  however,  screen/film  mammography  remains  a  diagnostic 
imaging  modality  where  image  interpretation  remains  very  difficult.  Breast  radiographs  are 
generally  examined  for  the  presence  of  malignant  masses  and  indirect  signs  of  malignancy 
such  as  the  presence  of  microcalcifications  and  skin  thickening.  Unfortunately,  it  is  unlikely 
that  major  improvements  in  imaging  performance  will  be  achieved  by  technical  advances  in 
screen/film  radiography  alone. 

The  major  reason  for  poor  visualization  of  small  malignant  masses  is  the  minor 
difference  in  x-ray  attenuation  between  normal  glandular  tissues  and  malignant  disease  [5]. 
This  fact  makes  the  detection  of  small  malignancies  problematical,  especially  in  younger 
women  who  have  denser  breast  tissue.  Although  calcifications  have  high  inherent 
attenuation  properties,  their  small  size  also  results  in  a  low  subject  contrast  [6].  As  a 
result,  the  visibility  of  small  tumors,  and  any  associated  microcalcifications,  will  always  be 
a  problem  in  mammography  as  it  is  currently  performed  using  analog  film. 

Improvements  in  the  ability  of  screen/film  mammography  to  detect  small  tumors  and 
microcalcifications  is  more  likely  to  occur  by  improving  the  visibility  of  these  features.  It 
has  been  suggested  that  as  normally  viewed,  mammograms  display  only  about  3  %  of  the 
information  they  detect!  [7]. 

Our  approach  to  feature  analysis  and  classification  is  motivated  in  part  by  recently 
discovered  biological  mechanisms  of  the  human  visual  system  [8].  Both  multiorientation 
and  multiresolution  are  known  features  of  the  human  visual  system.  There  exist  cortical 
neurons  which  respond  specifically  to  stimuli  within  certain  orientations  and  frequencies. 
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In  this  report  we  describe  exciting  new  results  accomplished  during  our  third  year  of 
study.  In  addition,  we  have  continued  our  efforts  in  the  development  of  wavelet  transforms 
that  exploit  orientation  and  frequency  selectivity  to  make  mammographic  features  more 
obvious  through  localized  contrast  gain.  Below,  we  describe  in  executive  summary,  recent 
accomplishments  related  to  the  objectives  stated  in  Phase  III  and  Phase  IV  of  our 
Statement  of  Work  (SOW).  The  methodology,  approach  and  experimental  studies  are 
reported  in  detail  in  the  body  of  this  report. 

1.1  Overview  of  Contents 

In  the  next  section,  we  describe  in  detail  our  wavelet  processing  algorithms,  experimental 
methods  and  example  results  obtained.  In  addition,  we  list  and  summarize  publications  of 
our  researchers  during  the  past  year  of  our  investigation. 

We  present  our  underlying  criteria  for  designing  a  database  of  105  abnormal  cases 
which  are  used  in  the  evaluation  of  our  computer-based  methods  for  enhancement, 
detection  and  analysis  of  masses.  We  take  into  consideration  a  number  of  issues  related  to 
image  acquisition,  image  selection,  and  image  annotation,  and  discuss  how  each  of  these 
points  affects  our  design  criteria. 

We  present  a  fast  implementation  of  the  discrete  dyadic  wavelet  transform,  which  we 
have  extensively  used  for  the  enhancement  of  digital  mammograms  [9,  10,  11,  12].  The 
algorithm  takes  advantage  of  symmetries  of  filter  coefficients  in  the  filter  bank 
implementation  of  the  transform,  and  of  the  symmetry  of  a  mirror  extended  input  signal. 
The  described  approach  is  not  limited  to  image  sizes  that  are  powers  of  two,  and  therefore 
allows  extraction  of  arbitrarily  sized  rectangular  regions  from  digital  mammograms.  We 
derive  a  second  derivative  analog  to  Mallat’s  “first  derivative  of  a  spline”  wavelet  [13]  and 
compute  filter  coefficients  for  a  discrete  filter  bank  implementation. 

In  our  earlier  report,  a  discrete  dyadic  wavelet  transform  was  applied  successfully  for 
enhancement  of  mammographic  features.  This  was  an  octave-based  transform  with  its 
scales  increasing  as  a  power  of  2,  and  had  limited  capability  for  charaterizing  signal  details 
at  scales  between  two  consecutive  levels,  including  band-limited  features,  such  as 
microcalcification  clusters,  spicular  lesions,  and  circular  (arterial)  calcifications.  To 
alleviate  this  problem,  and  to  more  reliably  identify  features  while  isolating  noise  through 
scale  space,  we  have  developed  a  sub-octave  wavelet  transform.  This  is  a  generalization  of  a 
dyadic  wavelet  transform  and  has  a  better  time-frequency  resolution  for  mammographic 
feature  analysis.  Our  experimental  results  demonstrate  advantages  over  traditional  dyadic 
wavelet  transforms  for  mammogram  enhancement.  We  also  describe  a  noise  reduction 
mechanism  to  avoid  amplifying  noise,  such  as  scatter  or  digitization  noise,  during  image 
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enhancement  processing. 

Linear  phase  filters  have  proven  advantageous  in  signal  and  image  processing  because 
they  limit  image  artifacts.  We  applied  complex  Daubechies  wavelets,  with  symmetric 
compact  support,  for  extracting  mammographic  features.  Also,  we  present  a  novel 
approach  for  denoising  using  complex  Daubechies  wavelets.  We  trace  the  evolution  of 
features  across  distinct  scales  by  their  wavelet  coefficients.  By  virtue  of  linear  phase  and 
dual  representations  of  scaling  functions,  we  show  how  noise  can  be  removed  while  features 
preserved. 

Most  image  processing  and  computer  vision  algorithms  assume  knowledge  of  the  proper 
spatial  scale  in  the  image.  We  describe  an  automatic  technique  that  adaptively  sets  the 
proper  spatial  scale  based  on  the  size  of  local  features  in  an  image.  Our  current 
implementation  automatically  determines  the  proper  scale  based  on  a  dynamic  system 
model  (motivated  by  analog  circuit  concepts)  that  continuously  searches  scale  space.  In  the 
discussion  below,  we  discretize  and  oversample  a  scale  space  representation  and  perform  an 
optimization  over  this  finite  set. 

We  show  that  regions  corresponding  to  masses  can  be  identified  reliably  through 
multiscale  features  of  a  continuous  wavelet  transform.  In  contrast  to  traditional  methods  of 
wavelet  analysis,  which  compute  transform  coefficients  at  dyadic  levels  of  scale,  we  show 
that  subtle  characteristics  of  mammographic  features  require  a  finer  parameterization  of 
scale  space. 

Next,  in  this  report,  we  show  the  derivation  of  a  scheme  to  calculate  the  wavelet 
transform  for  an  arbitrary  scale.  In  addition,  closed  contours  of  zero-crossing  features  are 
extracted  from  a  continuous  scale  frame  representation,  providing  precise  spatial 
localization  for  supporting  local  processing  such  as  segmentation  and  adaptive  contrast 
enhancement. 

Digital  radiographs  of  mammographic  features  (fibrils  and  masses)  within  an  RMI 
phantom  were  processed  for  demonstration.  Our  results  show  that  a  mass  and  fibril  known 
to  exist  in  the  physical  RMI  phantom  but  not  visible  via  traditional  image  processing 
techniques  (including  window  and  leveling,  histogram  equalization,  contrast  stretching, 
etc.)  are  made  clearly  visible  after  processing. 

We  present  results  of  local  enhancement  experiments  applying  hard  thresholding  of 
wavelet  coefficients.  The  images  used  were  synthetically  generated.  Care  was  taken  for 
realistic  modeling  of:  i)  noise,  ii)  Point  Spread  Function  (blurring  of  the  image),  and  iii) 
masses  that  may  exist  in  a  digital  mammogram.  Our  enhancement  techniques  were 
assessed  using  both  analytic  and  psychophysical  methods. 
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2.1  Database  Development 

2.1.1  Introduction 

In  this  section  we  describe  our  strategy  for  the  design  of  a  database  of  digitized 
mammograms.  The  database  can  take  a  number  of  forms,  from  a  simple  structured  file 
system  to  a  content-addressable  database  providing  intelligent  access  to  the  data.  In  order 
to  determine  the  type  of  database  to  construct  we  took  into  consideration  its  utility  in  the 
scope  of  the  project  and  the  availability  of  local  resources,  as  well  as  a  number  of  issues 
related  to  image  acquisition,  image  selection,  and  image  annotation.  Below,  we  discuss 
each  of  these  points  and  how  they  affected  our  design  criteria. 

2.1.2  Motivation 

The  first  step  in  designing  an  image  database  is  determining  its  intended  use,  as  this  has 
direct  impact  on  the  database  design  criteria.  Our  motivation  for  the  development  of  a 
database  of  digitized  mammograms  was  the  need  for  evaluating  our  computer-based 
methods  for  enhancement,  detection,  and  analysis  of  mammographic  abnormalities  on 
representative  cases. 

2.1.3  Image  Acquisition 

Good  quality  mammograms  are  essential  for  the  construction  of  a  database  of  digitized 
mammograms.  Kimme-Smith  [14]  suggests  that  as  a  minimum  requirement,  films  for  a 
database  of  digitized  mammograms  should  be  supplied  by  facilities  accredited  by  the 
American  College  of  Radiology  (ACR).  In  addition,  Rangayyan  et  al,  [15]  recommend  that 
films  should  come  from  institutions  such  as  teaching  hospitals  where  audits  of  all  aspects  of 
mammography  are  performed  on  a  regular  basis. 

In  addition  to  quality  mammograms,  selection  of  a  suitable  spatial  and  contrast 
resolution  is  fundamental  for  the  adequate  acquisition  of  digitized  mammograms.  The 
digitizer  has  to  be  capable  of  measuring  accurately  the  complete  range  of  optical  densities 
found  in  mammographic  images.  Nishikawa  [16]  recommends  that  the  number  of  bits  per 
pixel  should  be  chosen  so  that  image  noise  is  adequately  measured.  A  digitizer  measuring 
optical  densities  in  the  range  of  0.05-3.0  with  12  bits  of  quantization  is  generally  accepted 
as  a  suitable  choice  [17,  18,  16]. 

The  selection  of  acceptable  spatial  resolution  for  a  digitizer  depends  on  the  underlying 
task,  e.g.,  detection  of  microcalcifications  and  detection  of  masses.  It  is  widely  accepted 
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[17,  18,  16]  that  mammograms  with  microcalcifications  should  be  digitized  with  a  pixel  size 
of  50  microns  or  smaller  and  mammograms  with  masses  should  be  digitized  with  a  pixel 
size  of  100  microns  or  perhaps  larger. 

Other  general  recommendations  are  to  digitize  all  cases  on  the  same  digitizer  or  on 
digitizers  with  identical  or  correctable  characteristics  [16]  and  to  provide  calibration  data 
such  as  the  transfer  function  from  X-ray  exposure  to  pixel  values  [18]. 

In  our  database,  mammograms  were  digitized  with  a  Lumiscan  75  scanner  (Sunnyvale, 
California)  at  116-micron  pixel  resolution,  12-bit  gray  scale  resolution  and  optical  density 
range  of  0.05-3.00.  Cases  in  our  database  were  supplied  by  Shands  Hospital  at  the 
University  of  Florida,  an  ACR  accredited  and  FDA  certified  institution. 

2.1.4  Image  Selection  and  Annotation 

It  has  been  suggested  by  several  researchers  [19,  17,  16,  15]  that  a  database  which  is  to  be 
used  for  enhancement,  detection  and  analysis  of  mammographic  abnormalities  should 
contain  a  collection  of  images  in  which  a  wide  range  of  variation  of  normal,  benign  and 
malignant  conditions  are  adequately  represented.  In  addition,  a  database  of  digitized 
mammograms  should  provide  annotated  images  showing  the  location,  and  boundaries  of  all 
abnormalities  [19,  17,  18,  16,  15].  Furthermore,  cases  in  the  database  should  be 
unequivocally  categorized  as  normal  or  abnormal  [18,  14,  16].  Abnormal  cases,  benign  or 
malignant,  can  be  obtained  through  biopsy  [18,  14,  16].  Mammograms  regarded  as  normal 
can  be  obtained  by  selecting  cases  for  which  it  is  known  that  no  cancer  has  developed  for 
two  or  more  years  following  mammographic  screening  [18,  14,  16].  Image  annotations 
should  be  made  by  multiple  radiologists  experienced  in  reading  mammograms  and  a 
subjective  evaluation  of  the  detectability  of  the  lesion  should  be  provided  [19,  16].  It  has 
also  been  suggested  [19]  that  annotations  should  be  hierarchically  organized,  e.g., 
individual  microcalcifications,  clusters  of  microcalcifications,  and  the  overall  extent  of  the 
abnormality  should  all  be  recorded  and  related.  Also,  Astley  [19]  suggests  that  criteria  for 
classifying  images  within  the  database  be  established  and  standardized,  and  recommends 
that  radiographic  projection,  glandular  pattern  type,  diagnosis,  size,  location  and  subtlety 
of  abnormalities  be  included  as  part  of  such  criteria. 

Another  issue  of  importance  is  the  views  that  should  be  stored  in  the  database.  In 
standard  mammographic  screening,  a  case  consists  of  two  views  of  each  breast 
(craniocaudal  (CC)  and  mediolateral  oblique  (MLO)),  for  a  total  of  four  images.  There  is 
considerable  evidence  that  the  combination  of  the  craniocaudal  view  and  the  mediolateral 
oblique  projection  detect  more  cancers  than  a  single  view  [20].  In  fact,  two  views  per 
breast  are  the  only  accepted  norm  in  the  United  States.  Two  views  per  breast  allow  for  the 
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detection  of  cancer  through  right/left  asymmetry  and/or  the  confirmation  or  dismissal  of 
cancer  by  finding  its  presence  in  the  other  view.  Researchers  [17,  21,  16,  15]  have 
recognized  that  computer  algorithms  for  the  detection  of  mammographic  abnormalities 
must  use  more  than  a  single  view  and  that  a  practical  mammography  database  ought  to 
contain  at  least  the  four  views  of  a  standard  mammographic  screening  exam  and  perhaps 
previous  and  subsequent  examinations  of  the  same  patient. 

Consideration  must  also  be  given  to  the  distribution  of  normal,  benign,  and  malignant 
occurrences  within  the  database  [14,  16,  15].  Nishikawa  [16]  recommends  that  in  a 
mammography  database  the  number  of  normal  cases  should  approximately  be  the  same  as 
the  number  of  cases  for  each  lesion  type.  Also,  within  each  lesion  category,  there  should  be 
approximately  an  equal  number  of  benign  and  malignant  cases. 

As  a  starting  point  we  have  confined  the  lesion  type  to  masses  (mainly  due  to  the  pixel 
resolution  (116  microns)  of  our  scanner  and  constructed  a  database  of  105  abnormal  cases 
(420  images)  with  a  total  of  118  biopsy  proven  masses.  Following  the  recommendations 
outlined  above  we  have  established  criteria  for  including  cases  in  the  database  based  on  the 
American  College  of  Radiology  Breast  Imaging  Reporting  and  Data  System 
(ACR-BIRADS).  Cases  are  included  according  to  the  location,  shape,  margins,  density, 
size,  and  surrounding  glandular  pattern  of  the  underlying  abnormality,  and  the 
mammography  finding  (screening  assessment),  and  diagnosis  (biopsy),  see  Figure  1.  This 
information  allows  us  to  retrieve  cases  that  have  similar  characteristics  in  terms  of  location, 
lesion  appearance  and  diagnosis.  We  have  annotated  all  abnormalities  using  XMam  (see 
Figure  2)  with  the  assistance  of  a  radiologist  experienced  in  reading  mammograms  and 
included  a  subjective  evaluation  of  the  detectability  of  each  lesion  based  on  the  following 
scale: 

1.  Extremely  subtle:  Dense  or  heterogeneously  dense  tissue  with  the  margins  of  the 
mass  not  visible.  Asymmetry  or  architectural  distortion  may  be  present. 

2.  Very  subtle:  Dense  or  heterogeneously  dense  tissue  with  the  margins  of  the  mass 
partially  visible  on  one  view. 

3.  Subtle:  Heterogeneously  dense  with  the  margins  of  the  mass  partially  visible. 
Calcifications  may  or  may  not  be  present. 

4.  Relatively  obvious:  Heterogeneously  dense  with  all  margins  of  the  mass  visible  or 
clustered  microcalcifications  present  with  margins  of  the  mass  only  partially  visible. 

5.  Obvious:  Fatty  tissue  with  all  margins  of  the  mass  easily  visible  on  two  views. 
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Figure  1:  Case  characterization  form. 
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MLO  window 


Figure  2:  XMam  showing  a  malignant  tumor  with  spiculated  margins. 
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Using  XMaiti  capabilities,  annotations  were  organized  hierarchically  by  recording  individual 
masses  on  each  view  and  relating  instances  of  the  same  mass  across  views. 

Currently  our  database  consists  of  78  malignant  masses  (66.1  %)  and  40  benign  masses 
(33.9  %)  distributed  according  to  subtlety  as  follows;  28.8  %  obvious,  28  %  relatively 
obvious,  25.4  %  subtle,  14.4  %  very  subtle,  and  3.4  %  extremely  subtle.  During  the  current 
year  we  will  collect  more  cases  in  the  very  subtle  and  extremely  subtle  categories  as  well  as 
cases  of  normal  mammographic  occurrences. 

2.1.5  Summary 

We  have  described  some  underlying  criteria  to  be  followed  when  designing  a  database  of 
digitized  mammograms.  Respecting  these  criteria  we  constructed  a  database  of  105 
abnormal  cases  (420  images)  with  a  total  of  118  biopsy  proven  masses.  All  cases  were 
supplied  by  Shands  Hospital  at  the  University  of  Florida,  an  ACR  accredited  and  FDA 
certified  institution.  Mammograms  were  digitized  with  a  Lumiscan  75  scanner  (Sunnyvale, 
California)  at  116-micron  pixel  resolution,  12-bit  gray  scale  resolution  and  optical  density 
range  of  0.05-3.00.  All  abnormalities  were  annotated  hierarchically  by  a  radiologist  from 
Shands  Hospital  experienced  in  reading  mammograms  using  an  inhouse  tool  called  XMam 
described  previously  in  our  1994  annual  report.  Annotations  include  the  location,  shape, 
margins,  density,  size,  surrounding  glandular  pattern,  mammography  finding  and  diagnosis 
based  on  the  standardized  American  College  of  Radiology  Breast  Imaging  Reporting  and 
Data  System.  In  addition,  we  developed  a  subjective  evaluation  of  the  detectability  of  a 
lesion  and  evaluated  all  lesions  based  on  this  scale.  The  current  database  consists  of  78 
malignant  masses  (66.1  %)  and  40  benign  masses  (33.9  %)  distributed  according  to 
subtlety  as  follows:  28.8  %  obvious,  28  %  relatively  obvious,  25.4  %  subtle,  14.4  %  very 
subtle,  and  3.4  %  extremely  subtle. 

2.2  A  Discrete  Dyadic  Wavelet  Transform  Implementation 

2.2.1  Introduction 

Discrete  non-redundant  wavelet  transforms  have  been  successfully  applied  previously  in 
image  compression  applications  [22,  23].  However,  existing  problems  in  the  analysis  of 
medical  signals  and  images  motivate  the  use  of  redundant  wavelet  decompositions 
presented  in  this  report. 

The  discrete  dyadic  wavelet  transform  discussed  below  is  one  example  of  a  redundant 
representation.  As  originally  proposed,  the  wavelet  was  a  first  derivative  of  a  smoothing 
function,  and  was  used  as  a  multiscale  edge  detector  to  obtain  a  translation- in  variant 


15 


A  t 


representation  consisting  of  edges  [13].  A  reconstruction  algorithm  to  approximate  an 
original  signal  from  its  multiscale  edge  coefficients  alone  was  described  in  [13,  24]. 

In  comparison,  our  previous  work  [9,  10,  11,  12]  made  no  attempt  to  obtain  a 
parsimonious  representation  from  a  discrete  dyadic  wavelet  transform.  Rather,  the 
transform  was  intentionally  highly  redundant.  This  redundancy  was  exploited  for  image 
enhancement  by  modifying  the  transform  coefficients  and  reconstructing.  Here,  we  describe 
an  efficient  implementation  of  the  discrete  dyadic  wavelet  transform,  in  which  we  present  a 
tool  for  contrast  enhancement  in  digital  mammography. 

When  digital  filtering  of  a  finite-duration  discrete  signal  is  performed  via  circular 
convolution  the  filter  will  act  on  both  ends  of  the  signal  simultaneously.  This  may  lead  to 
artifacts  near  both  ends  of  the  result.  In  image  processing,  mirror  extension  of  an  input 
signal  is  a  very  popular  method  for  alleviating  such  boundary  effects.  Here,  we  present  a 
fast  filter  bank  implementation  of  the  discrete  dyadic  wavelet  transform  which  takes 
advantage  of  the  fact  that  the  input  signal  to  the  filter  bank  is  mirror  extended.  For 
clarity,  we  describe  both  one  and  two-dimensional  transforms.  An  extension  to  higher 
dimensions  is  given  in  [25]. 

2.2.2  One-Dimensional  Discrete  Dyadic  Wavelet  Transform 

A  discrete  wavelet  transform  is  obtained  from  a  continuous  representation  by  discretizing 
dilation  and  translation  parameters  such  that  the  resulting  set  of  wavelets  constitutes  a 
frame.  The  dilation  parameter  is  typically  discretized  by  fixed  dilation  steps  and  the 
translation  parameter  by  integer  multiples  of  a  fixed  step  [26].  Unfortunately,  the  resulting 
transform  is  variant  under  translations,  a  property  which  makes  it  less  attractive  for  the 
analysis  of  signals,  including  digital  mammography. 

Sampling  the  translation  parameter  with  the  same  sampling  period  as  the  input 
function  results  in  a  translation-invariant,  but  redundant  representation.  The  dyadic 
wavelet  transform  proposed  by  Mallat  and  Zhong  [13]  is  one  such  representation.  Let  us 
begin  with  a  brief  review  of  the  properties  of  a  dyadic  wavelet  transform  as  described  in 
[13],  included  here  for  completeness. 

The  dyadic  wavelet  transform  of  a  function  f{x)  G  L‘^(R)  is  defined  as  a  sequence  of 
functions 

where  W^fix)  =  f*  i^mix)  =  JZo  /(O  Mx-t)  dt,  ^p{x)  is  a  wavelet,  and 

■tfruix)  =  is  a  wavelet  expanded  by  a  dilation  parameter  (or  scale)  2"*.  Note 

the  use  of  convolution  instead  of  an  inner  product. 
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To  ensure  coverage  of  the  frequency  axis  the  requirement  on  the  Fourier  transform  of 
is  the  existence  of  >  0  and  Bi  <  oo  such  that 

OO 

A.  <  E  |^(2"w)p  <  B. 

m=— OO 

is  satisfied  almost  everywhere.  The  constraint  on  the  Fourier  transform  of  the  (nonunique) 
reconstructing  function  x(a;)  is  [13] 

OO 

5;  ^(2”c)x(2“u,)  =  l. 

m=— OO 

A  function  f[x)  can  then  be  completely  reconstructed  from  its  dyadic  wavelet  transform 
using  the  identity 

OO 

/(x)=  5:  W^f*Xm{x), 

m=— OO 

where  Xm{x)  =  2— x(2-™^). 

In  our  application,  processing  is  performed  on  discrete  rather  than  continuous 
functions.  When  the  function  to  be  transformed  is  in  the  discrete  form,  the  scale  2™  can  no 
longer  vary  over  all  m  G  a  finite  sampling  rate  constrains  the  window  size  from  being 
arbitrarily  small,  while  computational  resources  may  restrict  the  use  of  an  arbitrarily  large 
window.  Let  the  finest  scale  be  normalized  to  1  and  the  coarsest  scale  be  set  to  2^. 

The  smoothing  of  a  function  f{x)  €  L‘^{R)  is  defined  as 

Smf{x)  =  f  * 

where  =  2“’”(^(2~™x)  with  m  e  Z,  and  ^{x)  is  a  smoothing  function  (i.e.,  its 

integral  is  equal  to  1  and  <f>{x)  ^  0  as  |x|  — >  00  [13]). 

In  Mallat  et  al.  [13],  a  purely  real  smoothing  function  ^(x)  was  selected,  whose  Fourier 
transform  satisfied  ^ 

|^HP=  j:^(2’"a,)x(2”a-).  (1) 

m=l 

In  addition,  it  was  shown  that  any  discrete  function  of  finite  energy  {g{n)  €  1"^{Z))  can  be 
written  as  the  uniform  sampling  of  some  function  smoothed  at  scale  1,  i.e.,  g{n)  —  5'i/(n), 
where  /(x)  e  L'^{R)  is  not  unique.  Thus,  the  discrete  dyadic  wavelet  transform  of  Sif{n) 
for  any  coarse  scale  2^  is  defined  as  a  sequence  of  discrete  functions 

{SMfin  +  s),  {lT™/(n  +  s)}„e[i,M]}„g^, 

where  s  is  a  ^(x)  dependent  sampling  shift. 
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For  a  certain  choice  of  wavelets  the  discrete  dyadic  wavelet  transform  can  be 
implemented  within  a  fast  hierarchical  digital  filtering  scheme.  Next,  we  will  summarize 
the  relations  between  filters,  wavelets,  and  smoothing  functions. 

The  Fourier  transform  of  ^(x)  must  satisfy  [13] 

OO 

iiuj)  =  n  ^(2'^^),  (2) 

k=l 

where  j  stands  for  \/^,  the  lowpass  filter  frequency  response  H{uj)  is  differentiable,  and 

\H{u)\'^  +  \H{u>  +  7r)p  <  1  with  |if(0)|  =  1. 

Computing  Equation  1  for  the  finest  two  scales  shows  that 


0(2u;)  x{2u;)  =  \^{uj)\^-muj)\\ 

(3) 

while  computing  Equation  2  for  ^(2lo)  yields 

(4) 

If  we  choose 

(5) 

and 

x{2Lo)=e^^^K{uj)r{u), 

(6) 

where  and  are  digital  filter  frequency  responses,  *  denotes  complex 

conjugation,  and  insert  Equations  4,  5,  6  into  Equation  3  we  observe  a  relation  between  the 
filter  frequency  responses  [13], 

\H{uj)\^  +  G{Lo)K{io)  =  1. 

(7) 

Now  we  are  in  position  to  choose  filters  that  will  give  rise  to  wavelets  and  scaling 
functions  for  a  discrete  dyadic  wavelet  transform.  Filters  that  are  associated  with  a 
compactly  supported  orthonormal  wavelet  basis  are  certainly  a  possible  choice.  However, 
other  choices,  as  described  below,  provide  distinct  advantages  for  the  analysis  of  digital 
mammograms. 

Suppose  we  seek  a  wavelet  such  that  it  is  compactly  supported,  antisymmetric  or 
symmetric,  exhibits  good  edge  detection  capability  (i.e.,  equivalent  to  a  first  or  a  second 
derivative  of  some  smoothing  function),  and  is  as  regular  as  possible.  With  these  additional 
constraints,  orthonormal  wavelet  bases  are  almost  completely  ruled  out  (the  Haar  basis 
being  the  exception). 
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In  [10,  11],  Laine  et  al.  proposed  an  extension  to  the  family  of  filters  described  in  [13]. 
In  this  design,  the  wavelet  could  be  either  antisymmetric  and  equal  to  the  first  derivative  of 
some  smoothing  function  0(x),  or  symmetric  and  equal  to  the  second  derivative  of  0{x). 

When  the  wavelet  i(^{x)  is  antisymmetric  around  zero  (i.e.,  an  odd  function) 
is  an  odd  imaginary  function,  and  when  the  wavelet  il){x)  is  symmetric  around  zero  (i.e.,  an 
even  function)  is  an  even  real  function.  The  function  is  even  and  real 

in  both  cases. 

For  the  wavelet  ip{x)  to  be  the  first  (second)  derivative  of  some  smoothing  function 
0{x)^  '(piyj)  must  have  a  first-order  (second-order)  zero  at  u;  =  0,  and,  therefore,  G{oj)  must 
have  a  first-order  (second-order)  zero  at  c<;  =  0. 

However,  even  after  satisfying  all  of  the  above  constraints  there  remains  a  large  number 
of  possible  choices  for  H((jj).  In  [10],  the  class  of  filters  from  [13]  was  extended  by  choosing 

H{i^)  =  [cos  (8) 

where  2n-\-p  e  N,n  e  Z,  and  p  E  {0, 1}.  Once  H{u;)  is  chosen,  the  product  G{u))K{lo)  is 
constrained  by  Equation  7.  For  example,  choosing 

G{uj)  =  e-’tp  ^4;  sin  (9) 

determines 

A»  =  -l(e--G(c.))^  E  (cos(l))  .  (10) 

Note  that  H{lo)  is  a  lowpass  filter,  G{lo)  a  highpass  filter,  and  K{lo)  a  highpass  filter  for 
p  =  1  and  a  lowpass  filter  when  p  =  0. 

The  filters  described  in  Equations  8  through  10  are  finite  impulse  response  (FIR)  filters. 
For  implementation.  Table  1  and  Table  2  list  the  filter  coefficients  for  the  cases 
2n+p  e  {1,2}  and  2n+p  E  {3,4},  respectively. 


n  = 

=  0,  p  =  1 

n 

=  1 

,  p  =  0 

m 

h 

9 

k 

h 

9 

k 

-1 

0.5 

2 

0.25 

4 

-0.015625 

0 

0.5 

-2 

-0.125 

0.5 

-8 

-0.09375 

1 

0.125 

0.25 

4 

-0.015625 

Table  1:  Impulse  responses  of  filters  G{(jij),  and  K(u>)  for  2n+p  E  {1,2}. 

For  the  wavelet  to  be  symmetric/ antisymmetric  around  zero,  the  sampling  shifts  for 
this  family  of  filters  must  be  equal  to 
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n 

=  1 

,  p=  1 

n  — 

O 

II 

m 

h 

9 

k 

h 

9 

k 

-3 

-0.0009765625 

-2 

0.125 

-0.0078125 

0.0625 

-0.009765625 

-1 

0.375 

2 

-0.0546875 

0.25 

4 

-0.0458984375 

0 

0.375 

-2 

-0.171875 

0.375 

-8 

-0.13671875 

1 

0.125 

0.171875 

0.25 

4 

-0.0458984375 

2 

0.0546875 

0.0625 

-0.009765625 

3 

0.0078125 

-0.0009765625 

Table  2:  Impulse  responses  of  filters  G{u)),  and  K{lo)  for  2n+p  G  {3,4}. 


By  inserting  Equation  8  into  Equation  2  and  using  H^i  cos(2  ^o;)  =  ^  we  can  obtain 


(i>{uj)  = 


sin 


I  ) 


while  by  applying  Equation  9  and  Equation  5  we  see  that 

/  •  /  o;  \  \  2n-f-2 

2-. 


i,{^)  =  (jwf-”  ^  j  . 

Thus,  the  wavelet  'fpix)  is  a  first  (p  =  1)  or  a  second  (p  =  0)  derivative  of  a  smoothing 
function  0{x),  whose  Fourier  transform  is 

(*  /  u;  \  \  2n+2 

Note  that  d{x)  is  a  spline  of  degree  2n  +  1.  By  increasing  n,  e{uj)  becomes  more  localized  in 
the  frequency  domain  and  has  larger  support  in  the  spatial  domain. 

For  exposition,  Figure  3  shows  Oi^x^  for  n  G  {0,1,2}  and  the  corresponding  wavelets 
V>(a:)  for  2n+p  G  {1,2, 3, 4}.  With  the  exception  of  the  Haar  wavelet  (2n+p  =  1),  these 
wavelets  are  not  orthogonal  in  the  sense  that  their  binary  dilations  and  dyadic  translations 
generate  a  basis  of  L'^{R).  Their  regularity  order  is  2n+p-l  (i.e.,  ^(a;)  G 

Similar  to  orthogonal  and  biorthogonal  discrete  wavelet  transforms  [26],  the  discrete 
dyadic  wavelet  transform  can  be  implemented  within  a  hierarchical  filtering  scheme.  Using 
the  definition  of  the  discrete  dyadic  wavelet  transform  along  with  Equations  4  and  5  we  can 
formulate  the  analysis  section  of  such  a  filter  bank.  The  synthesis  section  simply  follows 
from  Equation  7.  Suppose 

FAu)  =  (11) 


where  F{u))  is  either  H{u\  G{u}\  or  K{lo')  (Equations  8  through  10).  We  may  then 
construct  a  filter  bank  implementation  of  the  discrete  dyadic  wavelet  transform  as  shown  in 
Figure  4. 
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Figure  3:  (a)  Primitives  0{x):  piecewise  linear  spline  (dashed),  cubic  spline  (solid),  and 
quintic  spline  (dotted);  (b)  wavelets  the  first  derivative  of  the  piecewise  linear 

spline  (dashed)  and  the  first  derivative  of  the  cubic  spline  (solid);  (c)  wavelets  xl^{x)  = 
the  second  derivative  of  the  cubic  spline  (dashed)  and  the  second  derivative  of  the  quintic 
spline  (solid). 
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Figure  4:  Filter  bank  implementation  of  a  one-dimensional  discrete  dyadic  wavelet  transform 
decomposition  (left)  and  reconstruction  (right)  for  three  levels  of  analysis.  denotes 

the  complex  conjugate  of 


Filters  referred  to  in  Equations  8  through  10  at  level  ?u-|-l  (i.e.,  filters  applied  at  some 
scale  2”^)  become  F(2’”a;),  where  F(u))  denotes  any  of  the  three  filters  at  level  1.  In  the 
spatial  domain  this  is  equivalent  to  upsampling  the  filter  impulse  response  by  2™  (i.e., 
inserting  2”^  -  1  zeros  between  subsequent  filter  coefficients  at  level  1).  Noninteger  shifts  at 
level  1  for  the  class  of  filters  with  p  =  1  are  rounded  to  the  nearest  integer.  An 
implementation  with  upsampling  of  filter  impulse  responses  (called  “algorithme  a  trous  ) 
was  first  proposed  by  Holschneider  tt  al.  [27].  The  complexity  of  such  a  filter  bank 
implementation  increases  linearly  with  the  number  of  levels. 

Let  us  refer  to  the  filters  at  level  1  (Equations  8  through  10)  as  “original  filters,”  to 
distinguish  them  from  their  upsampled  versions.  Let  an  input  signal  x[n)  be  real, 
x{n)  G  l\Z),  n  G  [0,iV  -  1],  and  let  X{uj)  be  its  Fourier  transform.  Depending  on  the 
length  of  each  filter  impulse  response,  filtering  an  input  signal  may  be  computed  either  by 
multiplying  X{oj)  by  a  filter  frequency  response  or  by  circularly  convolving  x{n)  with  a 
filter  impulse  response.  Of  course,  such  a  periodically  extended  signal  may  change  abruptly 
at  the  boundaries  causing  artifacts.  A  common  remedy  for  such  a  problem  is  realized  by 
constructing  a  mirror  extended  signal 


{  a;(— n  — 1)  if  n  G  [— A',  — 1], 
\  x{n)  if  n  G  [0,  A  -  1], 


(12) 


where  we  chose  the  signal  Xme(^)  to  be  supported  in  [— A,  A  —  1].  As  will  be  evident 
shortly,  mirror  extension  is  particularly  elegant  in  conjunction  with 
symmetric/antisymmetric  filters. 
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Let  us  first  classify  symmetric/antisymmetric  real  even-length  signals  into  four  types 
[28]: 

Type  I  /(n)  =  /(-n), 

Type  II  /(n)  =  /(-n  -  1), 

Type  III  fin)  -  -/(-n), 

Type  IV  fin)  =  -fi-n  -  1), 

where  n  €  [—N,N  —  1].  Note  that  for  Type  I  signals  the  values  at  /(O)  and  /(— iV)  are  not 
paired  with  any  other  value,  and  that  for  Type  III  signals  the  values  at  /(O)  and  fi-N) 
are  equal  to  zero. 

Using  properties  of  the  Fourier  transform  it  is  easy  to  show  that  the  convolution  of 
symmetric/antisymmetric  real  signals  results  in  a  symmetric/antisymmetric  real  signal.  If 
a  symmetric/antisynunetric  real  signal  has  an  even  length,  then  there  always  exists  an 
integer  shift  such  that  the  shifted  signal  belongs  to  one  of  the  above  types. 

Next  we  examine  the  filter  bank  implementation  of  a  one- dimensional  discrete  dyadic 
wavelet  transform  (Figure  4)  with  filters  derived  from  Equations  8  through  10  driven  by  a 
mirrored  signal  Xmei^)  at  the  input.  Let  the  number  of  levels  M  be  restricted  by 

4^-2 

M  <  log,  (13) 

-^max 

where  Rmax  is  the  length  of  the  longest  “original”  filter  impulse  response. 

Each  block  in  the  filter  bank  consists  of  a  filter  and  a  circular  shift  operator  (Equation 
11).  Equation  13  guarantees  that  the  length  of  the  filter  impulse  response  does  not  exceed 
the  length  of  the  signal  at  each  block. 

Since  an  input  signal  Xmein)  is  of  Type  II  and  noninteger  shifts  at  level  1  are  rounded 
to  the  nearest  integer,  it  follows  that  a  processed  signal  at  any  point  in  the  filter  bank 
belongs  to  one  of  the  types  defined  above.  This  means  that  filtering  a  signal  of  length  2N 
can  be  reduced  to  filtering  a  signal  of  approximately  one  half  of  its  length.  (For  Types  I 
and  III,  N  +  1  samples  are  needed.  However,  for  Type  III  one  needs  to  store  only  iV  -  1 
values  because  of  zero  values  always  stored  at  the  zeroth  and  (— V)-th  position  sample). 

Implementation  is  particularly  simple  for  filters  designed  with  p  =  0  (Equations  8 
through  10).  Filters  are  of  Type  I  in  this  case,  so  the  signal  at  any  point  of  the  filter  bank 
will  be  of  Type  II.  A  block  from  the  filter  bank  shown  in  Figure  4  (Equation  11)  can 
therefore  be  implemented  by 

R-1 

FMn)  =  fi0)unin)  +  ^  /(r)[u„(n  -  2™r)  +  unin  +  2™r)],  n  €  [0,  iV  -  1], 

r=l 
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where 


(14) 


Ull{l)  —  \ 


u{—l  —  1) 
u{l) 

u{2N-l-l) 


if  I  €  [— y?  “!]) 
iile[0,N-l], 
if  I  e  [N,  f  ], 


it(n)  is  an  input  signal  to  a  block,  m  +  1  is  a  level  (i.e.,  2™  a  scale),  f{n)  is  an  impulse 
response  of  some  original  filter,  R  is  the  length  of  the  filter,  and  N  is  the  length  of  the 
input  signal  x{n)  to  the  filter  bank. 

A  filter  bank  with  the  above  implementation  of  blocks  and  signal  x{n)  at  the  input 
yields  equivalent  results  as  circular  convolution  for  Xmei^i)  as  defined  by  Equation  12.  In 
addition  to  requiring  one  half  the  amount  of  memory,  the  computational  savings  over  a 
circular  convolution  implementation  of  blocks  are,  depending  on  the  original  filter  length, 
three  to  four  times  fewer  multiplications  and  one  half  as  many  additions. 

A  similar  approach  can  be  used  for  filters  with  p  ~  1.  However,  things  get  slightly  more 
complicated  in  this  case,  because  the  filters  are  not  of  the  same  type  and  the  signal 
components  within  the  filter  bank  are  of  distinct  type.  As  a  consequence,  an 
implementation  of  blocks  that  use  distinct  original  filters  will  not  be  similar,  and  the 
implementation  of  blocks  at  level  1  differs  from  the  implementation  of  blocks  at  other  levels 
of  analysis. 

The  decomposition  blocks  at  level  1  can  be  implemented  by  the  two  formulas 


Gsu{n)  =  g{0)[uij{n  —  1)  —  n  e  [1,  iV  —  1], 


Hsu{n)^  h{r)[unin  -  r  -  1)  +  un{n  +  r)],  nG  [0,7V], 

where  is  defined  by  Equation  14,  g{n)  and  h{n)  are  impulse  responses  of  the  filters 

computed  from  Equations  9  and  8,  respectively,  R  is  the  length  of  the  corresponding 
impulse  response,  and  N  is  the  length  of  the  input  signal  x[n)  to  the  filter  bank. 

Thus  decomposition  blocks  at  subsequent  levels  can  be  implemented  by 

Gsu{n)  =  g{0)[ui{n  -  2^s)  -  ui{n  +  2^s)],  n  G  [1,7V  -  1], 


and 


where 


Hsu{n)  =  y)  h{r)[ui{n  — 

r=0 


2™(r  +  5))  +  uj(n  +  2”^(r  +  s))], 


ui{l)  = 


! 


u{-l) 

u{l) 

u{2N  -  0 


if  I  G  [0,  iV], 
if /€  [iV  +  1,^], 


n  G  [0,  TV], 


(15) 


(16) 
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M(n)  is  an  input  signal  to  a  block,  m  +  1  identifies  a  level  of  analysis,  g{n)  and  h{n)  are 
impulse  responses  of  the  filters  from  Equations  9  and  8,  respectively,  R  is  the  length  of  the 
corresponding  impulse  response,  s  =  |  is  the  sampling  shift,  and  N  is  the  length  of  the 
input  signal  x{n)  to  the  filter  bank. 

The  outputs  from  blocks  Gs{2'^oj)  and  m  €  [0,M  -  1],  for  p  =  1  are  of  Types 

III  and  I,  respectively.  Next,  the  reconstruction  blocks  at  level  1  can  be  implemented  by 
the  two  formulas 

K-su{n)  —  ^  k{r)[uiii{n  —  r  +  1)  —  uui{n  +  r)],  n  €  [0,  iV  —  1], 

r=l 


and 


where 


R 

2 


Hlu{n)  =  ^  h{-r)[ui{n  -  r  +  1)  +  u/(n  +  r)],  n  G  [0,  iV  -  1], 


r=l 


\  —u{  —  l)  if/G[— y,— 1], 

0  if  /  =  0, 

uni{^)  =  {  w(0  if  1 G  [1,  iV  —  1],  (17) 

0  if  /  =  A^, 

-u(2iV-/)  if /G  [iV  +  1,^], 

ui{l)  is  as  defined  by  Equation  16,  u{n)  is  the  input  signal  to  each  block,  k{n)  and  h{n)  are 
impulse  responses  of  the  filters  from  Equations  10  and  8,  respectively,  R  is  the  length  of  the 
corresponding  impulse  response,  and  N  is  the  length  of  the  input  signal  to  the  filter  bank, 
x{n). 

Note  that  both  outputs  from  blocks  K-s{>^)  and  H*{uj)  for  p  =  1  are  of  Type  II. 

The  reconstruction  blocks  at  subsequent  levels  can  be  implemented  by 


K-su{n)  =  ^  k{r  +  l)[um(n  —  2’”(r  +  s))  —  uuj{n  +  2'”(r  +  s))],  n  G  [0,  W], 


r=0 


and 

H*u{n)  =  Hsu{n), 

where  u///(/)  is  given  by  Equation  17,  m  +  1  identifies  each  level,  k{n)  is  an  impulse 
response  of  the  filter  from  Equation  10  and  R  is  its  length,  5  =  |  is  the  sampling  shift,  N  is 
the  length  of  the  input  signal  to  the  filter  bank,  a:(n),  and  Hsu{n)  is  given  by  Equation  15. 
For  completeness,  we  mention  that  the  outputs  from  blocks  K-s{2'^lo)  and  i7*(2™a;), 
m  G  [1,  M  —  1],  for  p  =  1  are  of  Type  I. 

When  we  compare  the  above  implementation  of  blocks  for  p  =  1  to  circular  convolution 
driven  by  a  mirrored  signal  Xme{n)  at  the  input,  we  observe  that  approximately  twofold  less 
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memory,  four  times  fewer  multiplications  and  one  half  as  many  additions  are  required.  (For 
Type  I  signals  an  additional  sample  has  to  be  saved  because  two  values  are  without  a  pair) . 

The  implementation  presented  in  this  section  performs  all  operations  in  the  spatial 
domain.  However,  one  could  also  implement  the  structure  shown  in  Figure  4  with  the  input 
signal  Xme{n)  (Equation  12)  in  the  frequency  domain.  For  short  filter  impulse  responses, 
such  as  those  given  in  Tables  1  and  2,  the  spatial  implementation  described  in  this  section 
is  certainly  more  efficient.  For  long  filter  impulse  responses,  however,  filtering  is  faster  if 
implemented  in  the  frequency  domain.  Indeed,  it  makes  sense  to  keep  filtering  with  G(u;), 
which  has  never  more  than  three  nonzero  coefficients,  in  the  spatial  domain.  Additional 
details  on  alternative  implementation  strategies  can  be  found  in  [29]. 

2.2.3  Two-Dimensional  Discrete  Dyadic  Wavelet  Transform 

The  dyadic  wavelet  transform  of  a  two-dimensional  function  f2{xi,X2)  €  is  defined 

as  a  set  of  functions  [13] 

{^m/2(^i>^2),W2/2(ari,a:2)}„g 


where 

W^f2{xi,X2)  =  h  *  V’L(^i5^2)  =  /X!^/2(^i5t2)t/’m(^i~^i5®2— ^2)dtidt2  for  I  =  1,2  and 
V’m(^i5®2)  =  2-^'^4>‘{2-"^xi,2-'^X2)  are  wavelets  tp‘ixi,X2)  expanded  by  a  dilation 
parameter  2”^. 

To  ensure  coverage  of  the  frequency  space  there  must  exist  an  A2  >  0  and  B2  <  00  such 
that 

00  Z 

4I2  <  E  El0'(2“cc;i,2'"a;2)|"<  ^2 

771=  — 00  1=1 

is  satisfied  almost  everywhere.  If  (nonunique)  functions  *2);  X^(2?i, 2;2)  are  chosen 

such  that  their  Fourier  transforms  satisfy 

00  2 

^  ^^'(2-a;i,2-a;2)x'(2“u;i,2’"^2)  =  l, 

m=— 00  /=! 

the  function  f2{xi,X2)  may  be  reconstructed  from  its  dyadic  wavelet  transform  by 

00  2 

f2ixi,X2)=  E  E^m/2  *XL(®15^2), 

m=—oo  1=1 

where  —  2“^™x^(2“'” 2:^15  2'''”x2). 

However,  when  processing  discrete  functions  the  scale  2™  may  no  longer  vary  over  all 
m  e  Z.  Let  the  finest  scale  be  normalized  to  1  and  the  coarsest  scale  set  to  be  2^.  Let  us 


26 


introduce  a  real  smoothing  function  (l>2{xi,X2)  such  that  its  Fourier  transform  satisfies 

oo  2 

l*(wi,U2)P  =  Y,  (18) 

m—1  1=1 

Here,  as  in  one  dimension,  a  finite  energy  discrete  function  (5'2(«i5«2)  £  can  be 

written  as  the  uniform  sampling  of  some  function  smoothed  at  scale  1:  g2{ni,n2) 

=  S'i/2(ni,n2)5  where  f2{xi,X2)  G  is  not  unique,  and  Smf2{xi,X2)  =  f2*  (l>mixi,X2). 

Thus,  the  discrete  dyadic  wavelet  transform  of  5'i/2(ni,n2)  for  any  coarse  scale  2^  is 
defined  as 

{SMf2{ni  +  s,n2  +  s),  W^f2{ni  +  s,n2  +  s),  Wl/2(ni  +  s,n2  + 

where  s  is  a  wavelet  dependent  sampling  shift. 

To  implement  a  two-dimensional  discrete  dyadic  wavelet  transform  within  a  fast 
hierarchical  digital  filtering  scheme  the  wavelets  were  chosen  to  be  separable  products  of 
one-dimensional  functions: 

‘fp^ixi,X2)  =  i){xi)2<t){2x2), 

‘4’^{xi,X2)  -  i^{x2)  2(j){2xi).  (19) 

where  (f>{x)  and  '4>{x)  were  chosen  as  described  in  Section  2.2.2  (recall  that  the  Fourier 
transforms  ^(w)  and  ip{2(jo)  were  defined  by  Equations  2  and  5,  respectively). 

From  Equations  19,  5,  and  2  we  may  write 

^'(2a;i,  2u;2)  =  (20) 

where  G{uj)  is  the  frequency  response  of  a  digital  filter.  Choosing 

X'(2u;i,2a;2)  =  e^‘^^^K{ur)L2{u;2)H^S*M, 

x"(2cui,2w2)  =  (21) 

where  K{uj)  and  T2(w)  are  digital  filter  frequency  responses,  we  may  compute  Equation  18 
for  the  finest  two  scales  by 

^^'(2c<;i,2w2)  X^(2wi,2a;2)  =  |<^2(^i,‘<^2)|^  —  | ^2(2^1, 20)2) p.  (22) 

/=i 

Inserting  the  terms  defined  by  Equations  21,  20,  and  4  with  ^2{^ii^2)  =  <^(wi)^(^2)  into 
Equation  22  results  in 

/F(u;i)G(a;i)L2(^2)  +  K{u)2)G{u22)L2{loi)  +  \H[uji)\^\H(uj2)\^  —  !•  (23) 
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Equation  23  represents  a  relation  between  the  frequency  responses  of  the  digital  filters  used 
to  implement  a  two-dimensional  discrete  dyadic  wavelet  transform  and  is  a 
two-dimensional  analog  to  Equation  7. 

Solving  Equation  23  for  L2{oj)  by  substituting  K{lo)G{lo)  from  Equation  7  yields  the 
closed  formula 

i2(u.)  =  i(l  +  |ifHP)  m 

In  Table  3  we  provide  the  filter  coefficients  for  L2{lo)  from  Equation  24  for 
2n+p  e  {1,2,  3, 4}. 


h 

n  =  0, 

n  =  1, 

n  =  1, 

n  =  2, 

m 

p  =  1 

p  =  0 

p=  1 

p  =  0 

-4 

0.001953125 

-3 

0.0078125 

0.015625 

-2 

0.03125 

0.046875 

0.0546875 

-1 

0.125 

0.125 

0.1171875 

0.109375 

0 

0.75 

0.6875 

0.65625 

0.63671875 

1 

0.125 

0.125 

0.1171875 

0.109375 

2 

0.03125 

0.046875 

0.0546875 

3 

0.0078125 

0.015625 

4 

0.001953125 

Table  3:  Impulse  responses  of  filters  L2{io)  for  2n+p  G  {1,2, 3, 4}. 

As  in  the  one- dimensional  case,  a  two-dimensional  discrete  dyadic  wavelet  transform 
can  be  implemented  within  a  fast  hierarchical  filtering  scheme.  The  filter  bank 
implementation  follows  from  Equations  20,  4,  and  23,  and  is  shown  in  Figure  5. 

The  implementation  of  filtering  with  L2{co)  in  two  dimensions  (Table  3)  is  simple. 

Suppose  our  two-dimensional  structure  is  of  size  x  N2  and  its  mirror  extension 
handles  boundary  effects.  The  number  of  levels  remains  restricted  by  Equation  13,  where 
N  =  min{A^/ ;  /  G  [1,  2]}  and  Rmax  is  the  length  of  the  impulse  response  of 

Since  only  separable  filters  are  used  in  the  filter  bank  implementation,  we  can  simply 
apply  the  one-dimensional  implementation  as  presented  in  Section  2.2.2.  Equations  for 
implementing  blocks  6^5(0;),  and  at  distinct  levels  and  values  of  p 

are  then  straightforward. 

As  previously  mentioned  in  Section  2.2.2,  the  computational  savings  when  p  =  0  for 
this  implementation  over  circular  convolution  and  a  mirror  extended  signal  are  (depending 
on  the  original  filter  length)  three  to  four  times  fewer  multiplications  and  one  half  as  many 
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Figure  5:  Filter  bank  implementation  of  a  two-dimensional  discrete  dyadic  wavelet  transform 
decomposition  (left)  and  reconstruction  (right)  for  two  levels  of  analysis.  H*{uj)  denotes  the 
complex  conjugate  of  Hs{i^). 


additions.  Of  course,  for  longer  filters  a  frequency  domain  implementation  may  be  more 
efficient. 

2.2.4  Contrast  Enhancement  in  Digital  Mammography 

In  mammography,  early  detection  of  breast  cancer  relies  upon  the  ability  to  distinguish 
between  malignant  and  benign  mammographic  features.  The  detection  of  small 
malignancies  and  subtle  lesions  is  often  difficult.  Contrast  enhancement  can  make  more 
obvious  unseen  or  barely  seen  features  of  a  mammogram  without  requiring  additional 
radiation. 

Within  a  discrete  dyadic  wavelet  transform,  a  framework  for  contrast  enhancement  was 
achieved  by  applying  a  (possibly  nonlinear)  function  (referred  to  as  “enhancement 
function”)  to  wavelet  coefficients 

+  s,n2  +  s),  Wlf2{ni  +  s,n2  -F  ^)}me[i,M]} reconstructing  an 
enhanced  image  with  modified  coefficients  [10]. 

In  [9]  it  was  shown  that  multiscale  contrast  enhancement  techniques  based  on  three 
multiscale  representations:  1)  discrete  dyadic  wavelet  transform,  2)  (^-transform,  and  3) 
hexagonal  wavelet  transform,  can  outperform  traditional  contrast  enhancement  techniques 
such  as  histogram  equalization  and  unsharp  masking  [30]  when  applied  to  subtle  features  in 
digital  mammography.  Furthermore,  in  [12,  11,  10],  it  was  shown  that  unsharp  masking 
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with  a  Gaussian  lowpass  filter  can  be  formulated  as  a  special  case  of  contrast  enhancement 
via  a  discrete  dyadic  wavelet  transform. 

The  enhancement  function  was  chosen  such  that:  1)  low-contrast  areas  were  treated 
more  aggressively  than  existing  areas  of  high-contrast,  2)  edges  were  not  blurred,  3) 
monotonically  increasing,  and  4)  antisymmetric  (i.e.,  odd).  A  simple  function  that  satisfies 
the  above  conditions  is  [11] 


E{x) 


'  X  -  {K-  l)T  ifx<-T 
<  Kx  if  |a:|  <  T 

X  +  {K-  1)T  ifx  >  T. 


(25) 


Figure  6  shows  the  enhancement  function  from  Equation  25  for  parameter  values  K  20 
and  T  =  1. 


Figure  6:  The  enhancement  function  (Equation  25  with  K  =  20  and  T  =  1). 


Filters  with  p  =  0  (Equations  8  through  10)  were  found  to  be  more  suitable  for  contrast 
enhancement  than  wavelet  filters  designed  with  p  =  1.  It  has  been  shown  that  edges  change 
their  position  (introduce  artifacts)  when  the  wavelet  coefficients  computed  for  p  1  are 
modified  by  a  nonlinear  enhancement  function  [12]. 

Figure  7  shows  an  original  mammographic  image  and  an  enhanced  image  obtained  by 
using  filters  and  wavelets  with  2n-t-p  =  4  and  an  enhancement  function  defined  by  Equation 
25  with  K  =  30  and  T„  =  0.12max{|lF,(,/2(ni  +  3, ^2  +  5)| ;  I  €  {1,2},  (ni,n2)  G  Z^}, 
m  €  {1,  Af}.  The  image  matrix  size  was  775  x  436  pixels,  210  micron  spot  size,  and  the 
analysis  was  performed  up  to  the  fifth  level  (i.e.,  M  =  5). 

Enhancement  of  a  suspicious  area  of  a  mammogram  using  wavelets  with  2n-l-p  =  2  and 
an  enhancement  function  given  by  Equation  25  is  demonstrated  in  Figure  8.  The 
parameters  were:  K  =  20,  Tm  =  0.6max{|lF,(j/2(ui  +  -s,  +  -s)! ;  ^  G  {1,2},  (ui,  722)  G  ^  }, 

m  G  {1,M},  and  M  =  7.  The  image  and  enhanced  area  matrix  sizes  were  655  x  615  pixels 
and  175  x  165  pixels,  respectively. 
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Additional  details  on  contrast  enhancement  by  a  discrete  dyadic  wavelet  translorm, 
including  incorporation  of  denoising  into  the  enhancement  scheme,  can  be  found  in 
[10,  11,  12]. 


Figure  7:  (a)  An  original  mammographic  image  containing  a  spicular  mass,  (b)  An  enhanced 
image  with  spicular  borders  well  delineated. 
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2.2.5  Summary 


We  computed  the  coefficients  for  filters  of  low  order  splines  used  as  smoothing  functions  in 
a  fast  hierarchical  digital  filtering  implementation  of  a  discrete  dyadic  wavelet  transform 
for  both  one  and  two-dimensional  processing.  A  mirror  extended  input  signal  to  a  filter 
bank  implementation  of  the  discrete  dyadic  wavelet  transform  enabled  us  to  take  advantage 
of  the  symmetry/antisymmetry  of  both  signals  and  filters  for  an  efficient  implementation  of 
a  dyadic  transform.  There  are  no  restrictions  on  the  size  of  the  input  signal  which  makes 
this  scheme  attractive  for  extracting  arbitrary  rectangular  regions  from  digital 
mammograms.  The  presented  algorithm  is  iterative  across  scales:  a  noniterative  scheme 
suitable  for  high  speed  (interactive)  parallel  processing  can  be  derived  easily  by  using  the 
described  approach  on  each  dyadic  scale  independently. 

2.3  Denoising  and  Enhancement  via  Multiscale  Sub-Octave 
Analysis 

2.3.1  Introduction 

Orthonormal  wavelet  transforms  (OWT)  and  discrete  dyadic  wavelet  transforms  (DWT) 
have  been  successfully  used  in  various  application  areas,  such  as  data  compression,  edge 
detection,  texture  analysis,  noise  reduction,  and  image  enhancement  [31,  26,  32,  33,  13,  24]. 
In  this  part  of  our  research,  both  OWT  and  DWT  were  studied  for  denoising  and  contrast 
enhancement.  Analysis  and  experiments  show  that  a  DWT  has  several  advantages  over  an 
OWT  for  denoising  and  feature  enhancement  purposes  [34,  35].  A  DWT  with  a  wavelet 
equal  to  a  first  order  derivative  of  a  smoothing  function  as  its  basis  can  separate  feature 
variational  energy  (VE)  from  noise  VE  [34,  35].  This  helps  identify  which  wavelet 
coefficients  to  modify  to  enhance  specific  features  of  interest  (FOl)  without  amplifying 
noise.  The  mother  wavelet  was  a  smooth  and  antisymmetric  function  with  few  oscillations, 
which  kept  denoising  under  a  DWT  less  affected  from  pseudo-Gibbs  phenomenon.  This  can 
be  clearly  seen  from  the  results  under  an  OWT  [36,  37]  compared  with  the  results  obtained 
from  our  denoising  algorithm  under  a  DWT  [35].  The  discrete  filters  used  to  compute  a 
DWT  have  compact  support  requiring  few  taps.  In  addition,  DWT  wavelet  coefficients  have 
a  clear  meaning  in  that  they  are  proportional  to  the  image  intensity  changes  (gradients).  A 
DWT  is  a  stable  and  complete  representation,  which  guarantees  that  a  DWT  will  not  cause 
a  loss  or  gain  of  energy.  By  properly  modifying  signal  and  noise  VE,  we  can  achieve 
distinct  noise  reduction  and  feature  enhancement  within  a  digital  mammogram. 
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A  framework  for  multiscale  wavelet  analysis  was  implemented  as  a  research  tool  for 
denoising  and  feature  enhancement  [11,  12].  Here,  we  address  both  signal  or  image 
restoration  and  enhancement  problems  via  wavelet  shrinkage  and  feature  emphasis.  This 
approach  takes  advantage  of  both  soft-thresholding  and  hard-thresholding  wavelet 
shrinkage  techniques  to  eliminate  noise  while  preserving  the  sharpness  of  prominent 
features.  It  also  incorporates  local  energy  gain  for  nonlinear  processing  to  enhance  contrast 
within  structures  and  along  object  boundaries.  Thus,  feature  restoration  and  enhancement 
were  accomplished  by  modifying  the  gain  of  a  signal’s  VE.  Methods  have  been  developed 
for  achieving  denoising  and  feature  enhancement,  such  as  regulated  soft-thresholding  and 
adaptive  VE  gain  processing  combined  with  hard-thresholding  to  remove  small  noise 
perturbations  [34,  35]. 

Even  though  DWT-based  denoising  algorithms  have  produced  superior  results 
compared  to  other  methods,  we  point  out  that  the  DWT  has  limited  ability  to  characterize 
band-limited  features,  such  as  texture,  and  subtle  features  of  mammographic  images.  A 
DWT  is  an  octave-based  transform  with  its  scales  increasing  as  powers  of  2.  We  cannot 
visualize  or  detect  the  signal  details  for  scales  between  two  consecutive  levels  of  a  DWT.  To 
alleviate  this  problem,  and  to  more  reliably  isolate  noise  and  identify  features  through  scale 
space,  we  have  developed  a  sub-octave  wavelet  transform  [38] ,  which  is  the  generalization 
of  the  DWT.  A  sub-octave  wavelet  transform  provides  a  means  to  visualize  signal  details 
within  sub-octave  frequency  bands  with  equally- spaced  division  within  each  octave  band 
and  is  able  to  characterize  the  signal  details  more  precisely.  The  theoretical  development  of 
a  sub-octave  wavelet  transform  and  efficient  implementation  were  undertaken  in  this  study. 

2.3.2  Sub-Octave  Wavelet  Transform 

A  wavelet  transform  can  be  completely  characterized  by  a  family  of  wavelet  bases  used  for 
decomposition  and  reconstruction.  If  a  family  of  wavelets  bases  {i>n{x)}  is  complete  and 
orthonormal,  the  wavelet  transform  is  usually  called  an  orthonormal  wavelet  transform.  A 
Haar  wavelet  is  a  simple  example  of  orthonormal  wavelets.  The  problems  with  a  Haar 
wavelet  is  that  it  is  a  discontinuous  function  and  not  localized  in  the  frequency  domain. 
The  filters  {H,  G}  for  performing  an  orthonormal  wavelet  transform  satisfy  the  following 

equation 

\H{u)\^ +  \G{uj)\^  =  1.  (26) 

If  the  family  of  bases  is  complete  and  linearly  independent,  but  not 

orthonormal,  the  wavelet  transform  is  called  a  biorthogonal  wavelet  transform.  More 
generally,  if  the  family  of  wavelets  {i^n{x)}  is  not  linearly  independent  and  redundant,  they 

form  a  wavelet  frame  [39]. 
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For  a  discrete  dyadic  wavelet  transform,  we  relax  the  orthonormal  constraint  and  can 
have  different  decomposition  and  reconstruction  wavelets  as  long  as  the  corresponding 
high-pass  filters  G{ijj),K{u:i)  satisfy  the  following  equation 

\H{u)\^  +  G{u)K{u)  =  l.  (27) 

In  the  section  below,  we  develop  a  sub-octave  wavelet  transform  (SWT)  for  denoising  and 
enhancement. 


1-D  Sub-Octave  Wavelet  Transform 

In  this  section,  we  show  that  if  we  further  divide  each  octave  band  into  M 
equally-spaced  sub-octave  bands,  M  basis  wavelets  can  be  used  to  capture  the  detail 
information  of  a  signal  in  each  sub-octave  frequency  band.  The  M  wavelet  functions  are 
represented  as  €  L^{R),  where  m  -  1,2,  ..,M.  The  M-sub-octave  wavelet  transform 

of  a  1-D  function  f{x)  €  L‘^{R)  at  the  scale  2^'  and  location  x,  and  for  the  m-th  sub-octave 
frequency  band  is  defined  as 

W^f{x)  =  f*^^,{x),  (28) 

where  il>^,{x)  =  is  the  dilation  of  the  m-th  basis  wavelet  tl^'^{x),  sometimes  called 

the  m-th  analysis  wavelet,  at  scale  2\m-  1, 2, ...,  M,  and  j  €  Z.  In  the  frequency  domain, 
we  have 

=  /(cu)^™(2^'u;),  (29) 

by  taking  the  Fourier  transform  of  Equation  28.  To  have  a  more  flexible  choice  for  the  M 
basis  wavelets,  we  impose  that  the  wavelet  functions  satisfy  a  wavelet  frame  condition 
(similar  to  [13]) 

+00  M 

A  <  T.  H  <  B, 

j=—oo  m=l 

where  A  and  B  are  positive  constants  and  u;  G  i?.  In  the  spatial  domain,  we  have 

-f  CO  M 

AII/WII"  <  E  E  IIW'27/(^)II"  <  -B||/(*)IP- 


The  function  f{x)  can 
formula 


/(^)  = 


be  recovered  from  its  sub-octave  wavelet  transform  by  the  following 


+00  M  +00  M 

E  E  E/*«*72"W.  (30) 


j=—oo  m=l  j=-oo  m=l 

where  'y'^{x)  is  the  m-th  synthesis  wavelet  for  the  inverse  sub-octave  wavelet  transform. 
The  frequency  response  of  {t/>^(a;)  |  m  =  1, 2, ...,  M}  together  at  any  scale  2^  are  required  to 
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capture  the  details  within  an  octave  frequency  band.  For  perfect  reconstruction  of  f{x)^ 
and  7^(rc)  should  satisfy 

+00  M 

=  1.  (31) 

j=—Qo  m—1 

Equation  31  can  be  obtained  by  taking  the  Fourier  transform  on  both  sides  of  Equation  30. 
To  ensure  perfect  reconstruction,  the  whole  frequency  axis  should  be  covered  by  the 
analysis  and  synthesis  wavelets,  thus  the  reconstruction  wavelets  7™ (re)  can  be  any 
functions  whose  Fourier  transforms  satisfy  the  above  equation.  There  are,  in  fact,  many 
choices  for  the  decomposition  and  reconstruction  wavelets  to  satisfy  the  above  equation. 
Here  we  only  discuss  the  class  of  wavelets  which  are  the  first  order  derivatives  of  a 
smoothing  function,  such  as  spline  functions  of  any  order.  A  1-D  sub-octave  wavelet 
transform  can  be  easily  extended  to  a  2-D  sub-octave  transform  by  computing  sub-octave 
wavelet  coefficients  along  horizontal  and  vertical  directions,  as  explained  in  the  following 
section. 

2-D  Sub-Octave  Wavelet  Transform 

The  M-sub-octave  wavelet  transform  of  a  2-D  function  f{x,y)  €  at  the  scale  2^ 

and  location  (x,y),  and  for  the  m-th  sub-octave  frequency  band  is  defined  as 

Wlff{x,y)  =  f*4>irix,y),  (32) 

where  ^),  ^  =  1,2  (for  horizontal  or  vertical  direction), 

m  =  1, 2, ...,  Af,  and  j  G  Z.  In  the  Fourier  domain.  Equation  32  becomes 

The  function  f(x,y)  can  be  recovered  from  its  2-D  sub-octave  wavelet  transform  by  the 
following  formula 

+00  2  M 

fM=  Y.  EE 

j=—oo  m=l 
-foo  2  M 

=  E  EE  (3'‘) 

j=—oo  i=l  7n=l 

For  perfect  reconstruction  of  f{x,y),  i’^,{x,y)  and  'r^,{x,y)  should  satisfy 

+00  2  M 

E  EE  =  1.  (35) 

jz=—oo  i=l  m=l 

This  perfect  reconstruction  condition  can  be  obtained  by  taking  the  Fourier  transform  of 
Equation  34. 
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2.3.3  Discrete  Sub-Octave  Wavelet  Transform 

The  transform  parameters  in  a  sub-octave  wavelet  transform  are  continuous  variables.  This 
results  in  a  highly  redundant  representation.  It  is  possible  to  discretize  these  parameters 
and  still  be  able  to  achieve  perfect  reconstruction  [39].  For  digital  image  processing,  the 
sub-octave  wavelet  transform  of  a  discrete  function  can  be  carried  out  through  a  uniform 
sampling  of  a  continuous  sub-octave  wavelet  transform.  Below,  we  describe  the  discrete 
sub-octave  wavelet  transform. 

1-D  Discrete  Sub-Octave  Wavelet  Transform 

In  the  discrete  domain,  we  are  limited  by  the  finest  scale  of  1.  In  general,  a  discrete 
function  can  be  decomposed  into  the  details  of  fine-to-coarse  scales,  so  we  let 

-foo  M 

=  E  E  (i”(2'u,)7”(2'a.).  (36) 

m=\ 

For  a  J-level  discrete  sub-octave  wavelet  transform,  we  can  write 

j  M 

=  E  E  (2^'^)  +  (37) 

i=l  rn=\ 

If  the  scaling  approximation  of  a  function  f{x)  at  scale  2^  is  represented  as 

S2:if{x)  =  f*(P2j{^)^  (38) 

then  {W^f{xn),  S2jf{xn)  \  neZ,  j  =  1,2,...,  J,  and  m  =  1,2,...,M}  is  the  wavelet 
representation  of  a  discrete  function  /(a;„)  under  a  J-level  discrete  M-sub-octave  wavelet 
transform.  W^f{xn)  and  S2jf{xn)  are  uniform  samplings  of  their  continuous  counterparts 
respectively. 

If  we  let  J  =  1,  then  we  have 

M 

=  E  ^"‘(2a;)r(2u2)  +  |<^(2c^)p.  (39) 

m=l 

Let  us  assume  that  for  each  basis  wavelet  pair  ip‘^{Lo)  and  7'"(a;),  there  exists  a  pair  of 
corresponding  functions  and  /F™(u;)  which  need  to  be  determined  and  (with  certain 

temporal  shift  t)  satisfy 

i>^{2L0)  =  (^(tu)G™(u;)e-**‘^, 

7’"(2u;)  =(^(cu)/F'"(w)e*'^. 

For  scaling  function  there  exists  a  function  H{uj)  which  satisfies 

ip{2io)  = 
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Replacing  ^”"(2a;),  7™(2a;),  and  c^(2a;)  in  Equation  39,  we  have  the  following  fundamental 
relation  for  a  sub-octave  wavelet  transform  (SWT) 

M 

+  E  =  1.  (40) 

m—l 

The  discrete  sub-octave  wavelet  transform  of  a  function  f{xn)  €  can  be 

implemented  by  using  the  following  recursive  relations  between  two  consecutive  levels  j 
and  j  -f  1 

=  52,/(a;)G™(2M,  (41) 

4,+i/(u;)  =  S2Jicv)Hi2^uj),  (42) 

where  i  >  0,  1  <  m  <  M,  and  Sif{oj)  =  f{uj^,ujy).  The  reconstruction  5i/(w)  from  the 
sub-octave  wavelet  decomposition  can  be  implemented  through  the  following  recursive 
relation 

M 

Svf(^)  =  S2,*,f{oj)H’(2>u)  +  E  (43) 

m=l 

where  H*  is  the  complex  conjugate  of  H.  A  3-level  M-sub-octave  wavelet  decomposition 
and  reconstruction  process  based  on  the  above  recursive  relations  is  shown  in  Figure  9. 

The  corresponding  division  of  the  frequency  band  is  given  in  Figure  10.  In  general,  for 
M-sub-octave  decomposition  and  reconstruction  of  a  SWT,  we  need  to  have  M  pairs  of 
decomposition  and  reconstruction  basis  wavelets.  When  M  is  a  power  of  2,  we  can  perform 
the  decomposition  and  reconstruction  using  a  single  pair  of  wavelets  similar  to  the  way 
Daubechies  showed  in  [26].  Figure  11  presents  how  to  carry  out  a  2-level  4-sub-octave 
decomposition  and  reconstruction  of  a  SWT  using  FIR  filters  corresponding  to  a  single  pair 
of  basis  wavelets. 


2-D  Discrete  Sub-Octave  Wavelet  Transform 

For  the  decomposition  of  a  2-D  discrete  function,  we  let  the  frequency  response  of  a 
scaling  function  be  defined  in  the  following  formula 

+00  2  M  ^ 

|(^(a;i;,a>y)p  =  ^  E  ^  ‘x)y).  (44) 

j=l  i=l  m=l 

For  a  2-D  J-level  discrete  sub-octave  wavelet  transform,  we  can  formulate 


J  2  M 


|(,a(u;„u;,)|2  =  EE  E  ^‘’'"(2^u)„2H)f-'"(2^a>,,2^’cu,)  -h  |^(2-^l^„ 2^cc.,)i^  (45) 

j—1  izzl  m=l 

If  the  scaling  approximation  of  a  function  f{x,y)  at  scale  2-’  is  represented  by 

S23f{x,y)  =  f  *^2iix,y),  (46) 
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Figure  9:  A  3-level  SWT  decomposition  and  reconstruction  diagram  of  a  1-D  function. 


low  frequency  high  frequency 


Figure  10:  Division  of  the  frequency  band  under  the  SWT  shown  in  Figure  9. 


Figure  11;  A  2-level  4-sub-octave  decomposition  and  reconstruction  of  a  SWT. 


then  {  W^ffixu,  Vv),  S^jfixu,  Vv)  \u,v  e  =  1,2,  j  =  I, J,  and  m  =  1, M  }  is  called 
the  wavelet  representation  of  a  discrete  function  fi^XuiVv)  under  a  2-D  J-level  discrete 
M-sub-octave  wavelet  transform.  In  general,  (p{x,y)  is  a  2-D  scaling  function  and 
and  Y’”"{x,y)  are  the  m-th  directional  analysis  and  synthesis  wavelets.  There  are  many 
choices  of  these  functions  that  satisfy  Equation  45.  Similar  to  the  way  2-D  wavelets  are 
constructed  using  TD  wavelets  in  [13],  we  use  tensor  products  to  construct  2-D  sub-octave 
wavelets  using  TD  sub-octave  wavelets.  Thus  we  can  write 


^^’^{2^LO„2^iOy)  =  4,^(2^  LO,)cp{2^-^u;y), 

(47) 

(48) 

if{2^UJ^,2^LOy)  =  l^{2^LOx)‘P{2^0Jy). 

(49) 

Through  this  construction  of  2-D  wavelets,  we  implemented  a  2-D  sub-octave  wavelet 
transform  using  TD  convolution  with  FIR.  filters  of  a  TD  sub-octave  wavelet  transform. 
Figure  12  shows  the  division  of  the  frequency  plane  under  a  2-level  SWT.  A  2-D  sub-octave 
wavelet  transform  can  be  implemented  by  TD  FIR  filters  along  horizontal  and  vertical 
directions.  In  the  left  diagram  of  Figure  13,  a  fourth-order  spline  smoothing  function,  its 
first  and  second  derivative  sub-octave  wavelets,  and  the  corresponding  scaling  function  are 
shown.  The  right  diagram  of  Figure  13  presents  four  sub-octave  wavelets  corresponding  to 
the  first,  second,  and  third  derivatives  of  the  smoothing  spline  in  the  left  diagram. 
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A  dyadic  wavelet  transform  is  a  special  case  of  a  sub-octave  wavelet  transform  with 
M  =  1.  As  an  example,  a  discrete  2-sub-octave  wavelet  transform  is  shown  to  divide  the 
details  of  an  octave  band  into  the  details  of  2  sub-octave  bands,  and  one  sub-octave  can  be 
used  for  detecting  local  maxima  while  the  other  sub-octave  for  identifying  zero-crossings. 
Thus,  in  this  sub-octave  wavelet  representation,  we  take  advantage  of  local  maxima  and 
zero-crossing  representations  to  more  accurately  characterize  local  features. 

2.3.4  Noise  Modeling 

Signal  and  image  degradations  by  noise  and  artifacts  are  frequent  phenomena  of  the  real 
world  image  acquisition.  Image  degradations  have  a  significant  impact  on  radiologist 
performance  and  computer-assisted  methods.  In  addition,  noise  and  artifacts  often  make 
feature  extraction,  analysis,  and  recognition  algorithms  unreliable.  The  denoising  and 
feature  enhancement  techniques  described  next  improve  the  reliability  of  our  image 
processing  techniques. 

Additive  Noise  Model 

Additive  noise  can  be  (in  general)  represented  by  the  following  formula: 

/(x)  =  5f(x)  -F  7/„(x),  (50) 

where  5'(x)  is  an  unknown  function.  The  function  /(x)  is  a  noisy  observation  of  ^(x),  ?/a(x) 
is  additive  noise,  and  x  is  a  vector  of  spatial  locations  or  temporal  samples. 

For  additive  noise,  there  are  existing  techniques  based  on  mean  squared  error  or 
norm  optimization  to  remove  noise,  such  as  Donoho’s  wavelet  shrinkage  techniques 
[36,  40,  41],  Chen  and  Donoho’s  basis  pursuit  denoising  [42],  Mallat’s  local-maxima-based 
method  for  removing  white  noise  [24],  and  wavelet  packet-based  denoising  [43]. 

By  incorporating  a  denoising  mechanism  into  our  feature  enhancement  process,  we  may 
enhance  contrast  without  amplifying  noise.  Subtle  features  barely  seen  or  invisible  in  a 
mammogram,  such  as  microcalcification  clusters,  spicular  lesions,  and  circular  (arterial) 
calcifications,  may  be  enhanced  via  our  sub-octave  wavelet  transform.  Since  we  treat  noise 
and  features  independently,  we  were  able  to  obtain  superior  results  compared  to  existing 
algorithms  designed  for  enhancement  alone. 

2.3.5  Wavelet  Shrinkage  and  Feature  Emphasis 

In  this  section,  a  denoising  and  contrast  enhancement  scheme  based  on  wavelet  shrinkage 
and  feature  emphasis  is  presented.  In  this  investigation,  we  studied  hard  thresholding  and 
Donoho’s  soft  thresholding  wavelet  shrinkage  [40]  for  noise  reduction.  The  advantage  of 
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soft  thresholding  is  that  it  can  achieve  smoothness  while  hard  thresholding  better  preserves 
features.  In  our  unified  approach  for  accomplishing  denoising  and  feature  enhancement,  we 
took  advantage  of  both  thresholding  methods.  Donoho’s  soft  thresholding  method  [40]  was 
developed  for  an  orthonormal  wavelet  transform  [31]  primarily  with  Daubechies’s  Symmlet 
8  basis  wavelet  and  his  denoising  results  showed  some  undesired  side-effects,  from 
pseudo-Gibbs  phenomenon  [37].  By  using  a  SWT  and  basis  wavelets  without  any 
oscillations,  our  method  was  relatively  free  from  side-effects  after  denoising.  Experiments 
show  that  a  SWT  with  a  first  order  derivative  of  a  smoothing  function  as  its  basis  wavelet 
separated  feature  VE  and  noise  VE  coefficients.  In  this  algorithm,  we  adopt  regularized 
soft  thresholding  wavelet  shrinkage  to  remove  noise  VE  within  the  finer  levels  of  scale.  We 
then  apply  to  wavelet  coefficients  a  non-linear  gain  with  hard  thresholding  for  preserving 
features  while  removing  small  noise  perturbations  within  the  middle  levels  of 
space-frequency. 

Wavelet  Shrinkage  By  Soft  Thresholding 

Soft  thresholding  [40]  operation  can  be  represented  as 

u(x)  =  Ts(v(x),t)  =  sign(t;(x))  (|u(x)|  -  t)+,  (51) 

where  threshold  t  is  proportional  to  the  noise  level  and  x  G  the  domain  of  u(x),  and 
u(x)  is  the  result  of  soft  thresholding  and  has  the  same  sign  as  u(x)  if  non-zero.  The 
expression  (|z;(x)|  — 1)+  is  defined  as 

SWT  coefficients  can  be  processed  for  noise  reduction  by 

= Ts(w^rmAn,  (52) 

where  f  =  1, 2  (omitted  for  1-D  signals),  j  =  1, ...,  k,  k  <  J ,  m  —  1, ...,  M,  and  is  a  noise 
and  scale-related  threshold.  VE2}™’*/(x)  is  a  modified  coefficient.  We  recall  that  Donoho’s 
soft  thresholding  method  used  a  single  global  threshold.  However,  the  noise  coefficients 
under  a  SWT  have  average  decay  through  fine-to-coarse  scales,  we  regulated 
soft-thresholding  wavelet  shrinkage  by  applying  coefficient-dependent  thresholds  decreasing 
across  scales. 

Wavelet  Shrinkage  By  Hard-Thresholding 

A  hard  thresholding  operation  can  be  expressed  as 

u(x)  =  r//(u(x),f)  =  t;(x)(|u(x)|  >  t),  (53) 
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Soft  Thresholding  vs  Hard  Thresholding 


Figure  14:  Various  thresholding  methods. 


where  t  is  a  threshold  value,  x  £  D,  the  domain  of  'y(x),  and  u(x)  is  the  result  of  hard 
thresholding  and  has  the  same  sign  as  u(x)  if  non-zero.  The  meaning  of  (|t;(x)|  >  t)  is 
defined  by  the  expression 


(Kx)|>i)  = 


1  if  |u(x)|  >  <, 
0  otherwise. 


Sub-octave  wavelet  coefficients  were  modified  by  hard-thresholding  for  noise  reduction 

by 

ir‘r’7(x)  =  T„(wirf(x),  (‘f ),  (54) 

where  z  =  1,2  (omitted  for  1-D  signals),  j  =  and  k  <  J,  m  =  l,...,Af,  and 

(in  general)  a  threshold  related  to  noise  and  scale.  f(x)  is  a  processed  coefficient. 

Figure  14  shows  a  functional  diagram  of  normalized  soft  and  hard  thresholding. 


Feature  Emphasis  By  Generalized  Adaptive  Gain 

Recently,  several  wavelet-based  enhancement  techniques  have  been  developed 
[44,  45,  46,  9].  Adaptive  gain  non-linear  processing  [46,  10]  has  been  successfully  used  to 
enhance  features  in  digital  mammograms.  The  adaptive  gain  non-linear  processing  is 
further  generalized  in  this  study  to  incorporate  hard  thresholding  in  order  to  avoid 
amplifying  noise  and  actually  remove  small  noise  perturbations.  The  generalized  adaptive 
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gain  (GAG)  non-linear  operator  is  defined  as 


EoAoiv) 


'  0 

<  sign(n)  T2  +  a  (sigm(c(u  —  b))  —  sigm(— c(m  +  b))) 

V 


if  |t;|  <  Ti, 
if  T2  <  |n|  <  Ta, 
otherwise, 


(55) 


where  v  G  [—1,  1],  a  =  a(T3  —  T2),  u  =  sign(n)(|n|  —  T2)l{Tz  —  T2),  b  G  (0,  1), 
0  <  Ti  <  T2  <  Ta  <  1,  c  is  a  gain  factor,  and  a  can  be  computed  as 

_ 1 _ 

sigm(c(l  —  b))  —  sigm(— c(l  -|-  6))’ 


sigm(u) 


1 

1  -I-  e-^‘ 


(57) 


Ti,  T2,  and  T3  are  selected  parameters.  When  Ti  =  2^2  =  0  and  Ta  =  1,  The  expression  is 
equivalent  to  the  adaptive  gain  non-linear  function  in  [46,  10].  The  interval  [T2,  Ta]  serves 
as  a  sliding  window  for  feature  selectivity.  It  can  be  adjusted  to  emphasize  local  regions 
within  a  specific  range  of  variation.  Thus,  by  selecting  gain  values  and  local  windows,  we 
can  achieve  focused  enhancement  effects.  Through  this  non-linear  operator,  SWT 
coefficients  can  be  modified  for  feature  enhancement  by 


M'f  *  Egag 


(58) 


=  max 


(iwirfMi), 


(59) 


where  i  =  1,2,  1  <  m  <  M,  and  1  <  i  <  J.  W2)™’*/(a:,7/)  is  a  processed  coefficient.  Figure 
15  shows  a  functional  diagram  of  the  generalized  adaptive  gain  operator. 


2.3.6  Experimental  Results 

First  we  present  experimental  results  based  on  a  1-D  sub-octave  wavelet  transform.  This  is 
used  to  show  the  denoising  capability  of  our  method  with  features  restored  and  enhanced 
compared  with  other  well  known  denoising  methods.  Experimental  results  of  our 
SWT-based  denoising  with  features  restored  are  shown  in  Figure  16.  In  the  fourth  plot  of 
Figures  16(a)-(d),  the  denoised  results  are  overlaid  with  their  corresponding  originals. 

After  comparing  with  Coifman’s  and  Donoho’s  methods,  we  can  see  that  the  results  of  our 
SWT-based  method  are  visually  superior  and  relatively  free  from  side-effects  [37].  In  the 
next  experiment,  we  show  the  incapability  of  a  DWT  for  characterizing  band-limited  high 
frequency  features.  Figures  17(a)-(b)  show  the  original  and  noisy  “Doppler”  signals  with 
their  corresponding  5-level  DWT  and  a  coarse  scale  approximation.  The  finest  scale 
band-limited  features  (see  the  second  plot  in  Figure  17(a))  are  weak  and  obscured  when 
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Complex  Adaptive  Gain:  C=10,  B=0.35,  Tl=0.1,  T2=0.2,  T3=0.9 


Figure  15:  An  example  of  the  generalized  adaptive  gain  function. 

noise  is  present  (see  the  second  plot  in  Figure  17(b)).  These  high  frequency  band- limited 
features  are  lost  in  a  soft  thresholding  denoising  method,  which  is  shown  in  Figure  18(b). 
Figures  17(c)-(f)  show  a  2-sub-octave  wavelet  transform  of  the  original  and  noisy 
“Doppler”  signals.  Figures  17(c)  and  17(e)  show  first  sub-octave  coefficients  while  Figures 
17(d)  and  17(f)  display  second  sub-octave  coefficients.  Figure  18  shows  denoised  and 
enhanced  results  of  noisy  “Doppler”  under  a  DWT  and  a  SWT.  The  loss  of  high  frequency 
band- limited  features  and  other  factors,  make  the  result  from  the  DWT-based  method  less 
attractive  than  our  SWT-based  method. 

The  experimental  results  of  enhancement  with  noise  suppression  of  mammographic 
images  via  a  2-D  sub-octave  wavelet  transform  are  also  presented.  This  part  of  our 
experiment  demonstrates  the  potential  of  this  enhancement  method  without  amplifying 
noise,  such  as  background  noise  in  [46].  Figure  19(a)  shows  a  dense  mammographic  image. 
Enhancement  by  traditional  histogram  equalization  is  presented  in  Figure  19(b).  Figure 
19(c)  shows  the  result  of  SWT-based  enhancement  with  noise  suppression.  The  histogram 
equalization  method  makes  all  features  in  the  slightly  dark  areas  obscured.  The  mass  in 
the  middle  right  area  of  the  image  is  enhanced  and  becomes  visible  in  both  methods.  As 
shown  in  Figure  19(c),  enhancement  under  a  SWT  makes  barely  seen  tissue  structures 
visible.  Additional  enhancement  results  are  shown  in  Figure  20.  Figure  20  (b)  is  the  result 
presented  in  last  year’s  report  and  is  included  for  comparison.  Thus,  the  result  obtained 
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(a)  “Blocks 


(b)  “Bumps” 


Figure  16:  Denoised  and  restored  features  from  the  SWT-based  algorithm  from  top  to  bot¬ 
tom:  original  signal;  noisy  signal;  denoised  signal;  overlay  of  original  and  denoised  signal. 
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Figure  18:  Denoised  and  enhanced  results  of  noisy  “Doppler”  signal  under  a  DWT  and  a 
SWT  analysis. 


(c)  (d) 


Figure  20:  (a)  Mammogram  M56  with  blended  phantom  features,  (b)  Phantom  image, 
(c)  Nonlinear  enhancement  with  adaptive  wavelet  shrinkage  denoising  {Gm  =  20,  N  =  5, 
t  =  0.1).  (d)  SWT-based  enhancement  with  noise  suppression. 
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via  the  SWT-based  method  is  superior  to  that  of  traditional  DWT  methods  of  analysis. 

2.3.7  Summary 

We  have  shown  the  limitation  of  traditional  DWT  for  characterizing  band-limited  features, 
such  as  subtle  features  in  mammographic  images.  To  more  reliably  isolate  noise  and 
identify  features  through  scale  space,  we  formulated  and  developed  sub-octave  wavelet 
transforms.  If  the  number  of  sub-octave  bands  in  each  octave  is  a  power  of  2,  it  is  shown 
that  the  sub-octave  wavelet  transform  of  a  function  can  be  implemented  by  a  single  pair  of 
analysis  and  synthesis  wavelets.  FIR  filters  for  a  class  of  wavelets  designed  for  a  DWT  can 
be  used  to  implement  sub-octave  wavelet  analysis.  A  wavelet  shrinkage  scheme  for  noise 
suppression  was  presented  and  generalized  adaptive  gain  processing  was  developed  to 
enhance  features  without  amplifying  noise. 

Our  future  efforts  will  include  developing  feature  clustering  techniques  and  statistical 
analysis  for  efficiently  removing  noise  and  enhancing  features.  A  more  dedicated  feature 
cluster-based  technique  for  denoising  and  feature  enhancement  is  under  development.  We 
will  also  explore  iterative  methods  for  identifying  coherent  structures  from  an  image  via 
sub-octave  wavelet  analysis. 

2.4  Denoising  Using  Complex  Daubechies  Wavelets 

2.4.1  Introduction 

The  purpose  of  denoising  is  to  restore  an  image  that  differs  as  little  as  possible  from  its 
original.  Thus,  we  want  to  preserve  the  features  of  an  image  while  removing  noise. 

Here,  we  find  important  features  existing  at  one  scale  that  propagate  to  the  next  larger 
scale.  Wavelet  transforms  can  be  used  to  locate  the  singularities  of  an  image  because  of 
their  localization  characteristics.  By  taking  advantage  of  these  properties,  we  can  trace  the 
evolution  of  features  across  distinct  scales  by  their  wavelet  coefficients.  In  this  investigation 
we  present  a  new  approach  for  denoising  using  complex  Daubechies  wavelets.  By  virtue  of 
linear  phase  and  dual  representations  of  scaling  functions,  noise  can  be  removed  and 
features  preserved.  Preliminary  simulation  results  are  shown. 

2.4.2  Symmetric  Daubechies  Wavelets 

A  family  of  compactly  supported  orthonormal  wavelet  bases  were  first  constructed  by 
Daubechies  [31],  using  multiresolution  analysis.  In  oder  to  construct  a  compactly 
supported  wavelet  the  analysis  starts  from  a  scaling  function  <j)  with  compact  support. 

A  multiresolution  analysis  of  T^(R)  is  a  sequence  of  closed  subspaces  Vj  of  T^(R),  j  €  Z, 
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satisfying 

•  •  •  C  V2  C  Fi  C  Vo  C  V-i  C  VI2  C  •  •  • 

Since  the  scaling  function  0  G  Vo  C  V_i,  a  sequence  of  coefficients  hn  €  /^(Z)  exists  such 
that  the  scaling  function  satisfies 

<j){x)  =  V2YI  hn<f>{2x  -  n).  (60) 

n 

The  associated  wavelet  V’  can  be  generated  through  (j) 

i}){x)  =  ^'^gn<i>{2x  -  n),  (61) 


where 


Qn  =  {-inu 


with  h*  denoting  complex  conjugate  of  h. 

From  multiresolution  analysis,  f{x)  can  be  expressed  as 


Jo 


/(^)  ~  ^  '^Jo,n(^)  +  X/  ^ '*^nV’j>(^)5 

n  j—J  n 


(62) 


(63) 


where  jo  is  the  coarsest  level  scale  and  J  is  the  finest  level  scale.  The  coefficients  and  wi 
are  obtained  by  computing  the  inner  products 


4  =<  f’  ^i,k  >,  K=<  /,  i^j,k  >  ■ 


(64) 


Using  a  fast  wavelet  transform  (FWT)  algorithm,  we  can  compute  these  coefficients  by 
recursively  convolving  with  filters  h*^  and  g*,  respectively,  and  downsampling. 

<  =  ^12  K-2n4~^  >  <  =  V2  E  9l-2n4~'  >  (65) 

k  k 

where  is  the  original  signal.  The  original  signal  f{x)  can  be  reconstructed,  in  a  similar 
way,  via  the  inverse  FWT. 

4~^  =  V2'^hn-2ksi  +  V2'^9n-2kf^i-  (66) 

k  k 

The  block  diagrams  of  FWT  and  inverse  FWT  algorithms  are  shown  in  Figure  21. 
Multiresolution  analysis  leads  naturally  to  a  hierarchical  and  fast  scheme  for  the 
computation  of  the  wavelet  coefficients.  In  electrical  engineering  terms,  (65)  and  (66)  are 
the  analysis  and  synthesis  stages  of  a  subband  filtering  scheme  with  perfect  reconstruction. 
Filters  hn  and  gn  are  quadrature  mirror  filters  (QMF). 
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(b) 


Figure  21:  Block  diagram  of  fast  wavelet  transform:  (a)  decomposition;  (b)  reconstruction. 

A  linear  phase  property  can  be  achieved  with  symmetric  filters.  There  have  been  two 
distinct  methods  developed  to  generate  complex  Daubechies  wavelets.  Lawton  [47]  derived 
a  complex  valued  linear  phase  FIR  conjugate  quadrature  filters  (CQF)  from  a  class  of  real 
valued  FIR  CQF’s  by  replacing  certain  filter  zeros  by  their  reciprocal  conjugates.  Lina  et 
al.  [48]  used  two  parameterizations  to  compute  complex  solutions  for  Daubechies  wavelets. 
Let  the  Z-transform  of  hn  be 

H{z)  =  Y,Kz\  withF7(l)  =  l.  (67) 

n 

Symmetric  Daubechies  wavelets  have  the  following  advantages: 

1.  Compactness:  The  scaling  function  has  support  inside  interval  [-J,  J  +  1],  where  J  is 
the  order  of  the  underlying  analysis. 

2.  Orthogonality:  The  orthogonality  of  4>m,n  implies 

3.  Regularity:  The  derived  wavelets  will  optimize  the  regularity  of  a  given  compact 
support 

=  •  •  •  =  =  0. 

4.  Symmetry:  The  filter  must  satisfy  hn  =  h\-ni  i-e., 

H{z)  =  zH{z-^). 
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Figure  22:  Scaling  functions  and  wavelets:  (a),(b)  J=2;  (c),(d)  J=4 

The  standard  Daubechies  wavelets  have  the  first  three  properties,  while  the 
paxameterized  solutions  of  symmetric  Daubechies  wavelets  can  be  obtained  from  properties 
2,  3,  and  4.  Figure  22  shows  complex  scaling  functions  and  wavelets,  derived  by  Lina  et  al. 
[49,  48],  for  J  =  2  and  4,  where  (2J  +  2)  is  the  tap  number  for  the  complex  filters.  Note 
that  the  scaling  function  and  wavelet  for  J  =  4  are  smoother  than  those  for  J  =  2.  The 
real  part  and  the  imaginary  part  of  the  scaling  function  are  close  to  first  derivative  and 
second  derivative  of  a  Gaussian  function,  respectively. 

2.4.3  Methodology 
Previous  Techniques 

There  have  been  many  noise  reduction  methods  discussed  in  the  literature.  The 
simplest  noise  reduction  technique  is  an  equally  weighted  averaging  over  a  neighborhood 
(mean  filtering).  Averaging  produces  an  effect  similar  to  lowpass  filtering.  However, 
reducing  the  effects  of  random  noise  results  in  blurring.  The  method  of  median  filtering  is 
similar  to  mean  filtering  in  that  it  uses  the  local  median  instead  of  a  local  mean,  and 
exhibits  better  performance  for  edge-preservation  and  noise  reduction  (for  speckle  noise). 

More  recently,  denoising  based  on  wavelet  transforms  have  been  proposed.  Mallat  et  al. 
[24,  13]  used  the  local  maxima  of  the  wavelet  transform  modulus  to  analyze  image 
singularities.  The  algorithm  computed  Lipschitz  exponents  from  the  evolution  across  scales 
of  the  wavelet  transform  modulus  maxima,  removed  all  maxima  whose  amplitude  increased 
on  average  when  the  scale  decreased,  or  which  did  not  propagate  to  larger  scales.  A 
“denoised”  signal  was  then  reconstructed  with  an  alternative  projection  algorithm.  Mallat 
et  al.  [24,  13]  conjectured  that  a  signal  can  be  perfectly  reconstructed  from  the  locations 
and  values  of  the  modulus  maxima  of  the  wavelet  transform  at  different  scales.  This 
conjecture  has  been  disproved  by  Meyer  [50]  although  the  reconstruction  gives  a  good 
approximation  in  many  cases. 

Donoho  and  Johnstone  [40,  36]  proposed  a  three-step  method  for  recovery  of  a  noisy 
signal.  Their  method  is  based  on  a  nonlinear  shrinkage  of  wavelet  coefficients.  In  contrast 
to  the  hard  thresholding  of  wavelet  coefficients,  their  soft-thresholding  approach  provides  a 
smoother  means  for  removing  noise. 

Mallat  and  Zhang  [51]  introduced  a  matching  pursuit  algorithm  that  decomposed  a 
signal  into  a  linear  expansion  of  waveforms  that  belonged  to  a  redundant  dictionary  of 
functions.  With  a  dictionary  of  Gabor  functions,  a  matching  pursuit  defines  an  adaptive 
time-frequency  transform.  A  coherent  structure  is  chosen  by  testing  the  amplitude  ratio  of 
successively  extracted  amplitudes  with  thresholded  noise  components.  The  denoised  signal 
can  then  be  reconstructed  from  coherent  components.  There  is  a  trade-off  between 
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computation  time  and  precision.  If  we  use  too  many  bases  in  a  dictionary,  the  algorithm 
will  take  a  long  time  to  converge;  if  a  small  set  of  bases  is  used,  the  matching  may  not  be 
good  enough. 

Two  variations  of  matching  pursuit  algorithm  have  been  proposed.  Both  of  them 
express  a  signal  as  a  superposition  of  a  collection  of  libraries  of  waveforms.  The  first  one 
was  an  adapted  waveform  transform  by  Coifman  and  Wickerhauser  [52].  They  used 
wavelets,  wavelet  packets,  and  windowed  trigonometric  waveforms  as  the  libraries  of 
waveforms  because  these  bases  could  be  implemented  with  fast  algorithms,  and  thus  could 
achieve  high-speed  interactive  performance.  They  found  the  best  basis  in  each  library  by 
minimizing  a  cost  function.  A  signal  can  then  be  divided  into  two  parts  -  coherent 
components  and  incoherent  components  -  via  a  threshold.  The  incoherent  component  is 
then  considered  a  new  signal  and  decomposed  again.  This  process  is  performed  iteratively 
until  it  reaches  a  fixed  number  of  decompositions  or  no  coherent  components  are  obtained. 
The  denoised  signal  is  obtained  by  superimposing  the  coherent  components  together. 
Because  the  bases  are  orthogonal,  this  method  may  not  obtain  sparse  representations  if  a 
signal  is  composed  of  several  non-orthogonal  components. 

In  another  method,  basis  pursuit  was  presented  by  Chen  and  Donoho  [42].  They  used  a 
richer  libraries  of  waveforms  as  basis,  including  wavelets,  steerable  wavelets,  Gabor 
dictionaries,  multiscaJe  Gabor  dictionaries,  wavelet  packets,  cosine  packets,  and  chirplets. 
They  choose  the  basis  whose  coefficients  have  a  minimum  P  norm.  The  solution  can  be 
obtained  by  solving  an  equivalent  linear  program.  In  order  to  deal  with  a  noisy  signal,  they 
proposed  an  approximate  decomposition,  as  mentioned  previously,  by  solving  a  quadratic 
program.  Noise  in  both  variations  is  assumed  to  be  white  Gaussian. 

Freeman  and  Adelson  [53]  proposed  the  concept  of  steerable  filters  and  applied  it  to 
several  problems.  With  a  set  of  “basis  filters”,  one  can  adaptively  steer  a  filter  along  any 
orientation.  Hilbert  transform  pairs  were  constructed  to  find  a  local  “oriented  energy” 
measure.  Noise  was  then  reduced  by  performing  a  soft  threshold  function  [54]  on  a  pyramid 
of  multiscale  coefficients. 

It  is  well  known  that  linear  phase  filters  are  important  in  image  coding  [55].  There 
exists  no  real  orthogonal  compactly  supported  wavelet  bases  in  which  a  scaling  function  (f) 
has  linear  phase,  except  for  the  Haar  basis  [31,  39].  Complex  Daubechies  wavelets  were 
first  developed  by  Lawton  [47]  and  Lina  and  Marrand  [49,  48]  to  achieve  the  linear  phase 
property.  Lina  et  al.  proposed  a  sharping  algorithm  for  image  enhancement  based  on 
complex  Daubechies  wavelets  analysis. 
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Figure  23:  Block  diagram  of  ^'algorithme  a  trous’':  (a)  decomposition;  (b)  reconstruction. 


Complex  Daubechies  Wavelets  Approaches 

In  order  to  obtain  linear  phase,  we  can  relax  the  orthogonality  requirement  and  use 
biorthogonal  bases,  or  give  up  compactly  supported  bases  to  get  infinite  impulse  response 
(HR)  filters.  Also,  we  may  choose  to  drop  the  real  valued  filter  requirement  and  obtain 
complex  valued  linear  phase  FIR  filters.  We  found  that  complex  Daubechies  wavelets  have 
some  interesting  properties  which  may  be  useful  for  the  analysis  of  mammographic 
features,  such  as  dual  representation  -  local  extrema  and  zero- crossings.  We  have  used 
these  characteristics  in  our  denoising  algorithm. 

As  we  mentioned  previously,  for  compactly  supported  orthonormal  wavelet  bases,  we 
may  decompose  and  reconstruct  an  image  via  a  fast  wavelet  transform  algorithm  [56] . 
However,  this  transform  is  not  shift-invariant,  which  is  a  problem  in  mammographic 
analysis.  In  order  to  accomplish  shift-invariance,  we  needed  to  find  an  equivalent  structure 
without  downsampling.  We  used  the  structure  shown  in  Figure  23  to  implement  the 
transformation.  Compared  to  Figure  21,  every  filter  /  and  downsampling  (upsampling) 
operation  is  replaced  by  a  new  filter  D®/.  D®/  is  obtained  by  inserting  2®  —  1  zeros  between 
every  pair  of  the  coefficients  of  /.  That  is. 


P7I 


fn.  if  n  =  2®m, 

0  otherwise. 


where  m  and  n  are  lengths  of  the  original  filter  and  filters  at  z-th  level,  respectively.  This 
approach  is  called  ^^algorithme  a  trous”  (algorithm  with  holes)  [27,  57].  The 
decomposition  and  reconstruction  equations  at  level  i  are  as  follows. 


58 


Test  signal  R®al  components 

(a)  (b) 

Figure  24:  (a)  Test  signal,  (b)  Four  levels  of  waveform  basis. 

Decomposition: 

5*+^  =  {D^h*)*s\  (68) 

^i+i  ^  {D^g*)*s\  (69) 

Reconstruction: 

5*  =  {D^h)  *  +  {D^g)  *  (70) 

where  indicates  discrete  convolution. 

One-dimensional  case 

Edges  are  important  for  detecting  mammographic  features.  In  this  study,  we  considered 
two  kinds  of  features:  steps  (including  low  to  high  and  high  to  low)  and  pulses,  which  are 
shown  in  Figure  24(a).  Real  components  of  four  levels  of  wavelets  coefficients  were  also 
shown  in  Figure  24(b).  Basically,  each  level  had  two  types  of  waveforms.  (Actually,  there 
are  four  kinds  of  waveforms.  Two  waveforms  form  a  pair  and  each  pair  is  opposite  of  each 
other.)  The  waveforms  among  distinct  levels  have  similar  shapes,  except  for  dilations.  We 
saved  the  waveforms  in  the  fourth  level  and  normalized  them  as  a  basis.  Since  these 
waveforms  are  symmetric  or  antisymmetric,  we  only  needed  to  save  half  of  the  samples. 

We  can  decompose  a  signal  into  several  levels  of  wavelet  coefficients  using  complex 
Daubechies  wavelets.  The  filter  coefficients  that  we  used  are  shown  in  Table  2.4.3  for 
J  =  6.  Because  the  filters  are  complex,  the  coefficients  in  each  level  contain  two  parts:  real 
and  imaginary  components.  According  to  our  experiments,  four  levels  of  wavelet 
coefficients  were  enough  for  analysis.  If  we  decomposed  further,  some  of  the  features  were 
smoothed  and  disappeared. 

We  started  the  analysis  from  the  third  level.  Real  components  were  considered  first. 
The  coefficients  were  thresholded  and  local  extrema  were  identified.  The  features  of  a 
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Figure  25:  Four  levels  of  wavelet  coefficients  for  a  noisy  input  sample:  (a)  real  components; 
(b)  imaginary  components. 


(a)  (b) 

Figure  26:  Four  levels  of  denoised  wavelet  coefficients:  (a)  real  components;  (b)  imaginary 
components. 

signal  should  propagate  across  scales.  Extrema  in  the  third  level  were  compared  with  those 
in  the  fourth  level.  We  kept  the  extrema  that  also  appeared  at  the  fourth  level.  Extrema 
that  did  not  show  up  at  the  fourth  level  were  considered  as  false  peaks  and  removed. 

Because  the  lengths  of  the  basis  were  long,  four  simple  3-tap  masks,  {-1.0,  3.0,  -1.0}, 
{-0.1,  2.5,  -2.5},  {1.0,  -3.0,  1.0},  and  {0.1,  -2.5,  2.5},  were  used  to  classify  the  waveforms. 
These  four  masks  are  the  simplified  forms  of  the  basis  shown  in  Figure  24(b)  and  were  used 
to  identify  the  waveform  type.  At  the  third  level,  we  simply  aligned  the  center  of  the  mask 
with  the  extrema  and  the  other  two  of  the  masks  with  the  eighth  neighbors  of  the  extreme. 
These  were  then  multiplied  and  summed.  The  mask  which  had  a  maximum  value  was  used 
to  decide  the  waveform  type  of  a  feature. 

Since  we  classified  the  waveform  type,  we  traced  backwards  to  the  first  two  levels  and 
find  the  positions  of  extrema.  A  feature  will  propagate  across  different  scales  with  the  same 
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dc  component  dc  component 


Figure  27:  (a)  Noisy  dc  components,  (b)  Denoised  dc  components. 


waveform  type.  We  kept  the  best  match  at  the  first  two  levels.  Then  we  replaced  real 
wavelet  coefficients  in  four  levels  with  corresponding  idealized  waveform  basis  proportional 
to  the  magnitudes  of  the  extrema.  Because  the  basis  were  fitted  for  the  fourth  level,  we 
downsample  by  two  to  obtain  the  basis  for  the  other  levels.  Imaginary  components  of  the 
corresponding  periods  were  preserved.  The  reasons  are  as  follows:  the  energy  in  the 
imaginary  part  of  the  coefficients  is  significantly  smaller  than  that  in  the  real  part,  and 
coefficients  of  a  noisy  signal  are  similar  to  those  without  noise.  Four  levels  of  wavelet 
coefficients  for  noisy  input  are  shown  in  Figure  25,  and  the  denoised  versions  are  shown  in 
Figure  26. 

Next,  we  focus  on  the  DC  component.  This  component  is  obtained  from  passing  an 
input  through  several  lowpass  filters,  and  is  very  smooth.  We  already  know  that  the  real 
part  and  imaginary  part  of  a  scaling  function  is  close  to  the  first  derivative  and  the  second 
derivative  of  Gaussian,  respectively.  We  first  locate  regions  within  imaginary  components 
that  are  less  than  some  threshold  and  oscillate  near  zero.  We  set  these  regions  to  zero  and 
replace  each  corresponding  region  of  real  components  with  its  average  value.  Figure  27 
shows  a  sample  noisy  dc  component  and  denoised  version.  Finally,  we  reconstructed  the 
original  signal  with  the  modified  coefficients. 

The  algorithm  is  summarized  as  follows: 

1.  Using  a  frame  representation,  we  decompose  a  signal  into  several  levels  of  complex 
wavelet  coefficients. 

2.  Starting  from  the  third  level,  we  identify  the  local  maxima  for  real  wavelet 
coefficients  that  are  greater  than  some  threshold. 

3.  Selected  coefficients  are  compared  with  the  fourth  level  to  remove  unwanted  maxima 
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and  classify  each  waveform  type. 

4.  Trace  backward  to  the  first  two  levels  to  find  the  positions  of  maxima. 

5.  Replace  real  wavelet  coefficients  with  corresponding  idealized  waveform  basis  and 
preserve  imaginary  components. 

6.  Locate  regions  within  the  imaginary  component  of  DC  channel  that  are  less  than 
some  threshold  and  oscillate  near  zero. 

7.  The  real  part  of  each  corresponding  region  is  replaced  with  its  average  value. 

8.  Reconstruct  original  signal  using  modified  coefficients. 

Two-dimensional  case 

The  two-dimensional  case  is  easily  extended  from  the  one- dimensional  case  by  tensor 
product  construction.  In  order  to  reduce  computational  cost,  we  used  a  method  proposed 
by  Mallat  et  al.  [13]  to  find  features  along  the  x-axis  and  y-axis,  respectively.  Let  the 
scaling  function  ^{x)  be  expressed  by 

-1-00 

^{lo)  =  n  (71) 

p-1 

where  H{u)  is  the  Fourier  transform  of  a  lowpass  filter  h{n)  and  s  is  the  sampling  shift. 
This  implies  that 

^{2io)  =  (72) 

A  wavelet  rp{x)  whose  Fourier  transform  ^(cu)  can  be  defined  in  the  same  way: 

i>(2uj)  =  (73) 

where  G{u>)  is  the  Fourier  transform  of  a  highpass  filter  g{n).  The  two  wavelets  ij;^{x,y) 
and  'tp'^{x,y)  are  separable,  and  are  defined  as 

y)  =  2il){x)(f>{2y),  (74) 

'ip‘^{x,y)  =  2(j){2x)(j){y).  (75) 

The  Fourier  transform  of  these  two  wavelets  are  then  given  by 

'tp‘^{2uJx,2LOy)  =  ^{u!x)G{iL!y)^{LOy). 
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(76) 

(77) 


The  two-dimensional  wavelet  transform  of  f{x,y)  can  be  expressed  as  the  set  of  functions 


Wf  =  {W^J{x,y),W^J{x,y)), 


(78) 


where  j  €  Z.  We  now  find  two  reconstructing  functions  and  such  that 

f{x,y)  can  be  reconstructed  from  its  wavelet  transform 

+  00 

f{x,y)  =  (^2^/  *  xh{x, y)  +  *  xv{x,y))-  (79) 

j=-oo 

In  order  to  obtain  perfect  reconstruction,  the  Fourier  transforms  of  x^i^'>y)  and  x^{x,y) 
must  satisfy 

-f-CO  ^  ^ 

2^Uy)  +  2^ujy)x^{2^oj^.  2^'u;,))  =  1.  (80) 

j=-oo 

Let  and  x'^i^x^^y)  be  defined  such  that 

Xi(2w.,  2ujy)  =  (81) 

X2(2u;,,  2uy)  =  e^^‘^yK{cOy)LM^M^{ujy),  (82) 


where  K{u})  and  L{uj)  are  27r  periodic  functions  and  they  have  to  satisfy 


Giu)K{io)  +  \H{io)\'^  =  1, 


(83) 

(84) 


In  our  algorithm,  /(n)  and  g{n)  are  lowpass  and  highpass  filters,  respectively.  We 
choose  k{n)  equals  to  g*{n)  to  meet  the  first  criterion.  Next  we  need  to  generate  a  filter 
l{n)  which  will  satisfy  the  second  criterion.  In  Mallat’s  original  derivation,  he  developed  his 
algorithm  based  on  a  continuous  wavelet  transform.  We  approximate  his  derivation  with  a 
discrete  formulation. 

The  support  of  the  lowpass  filter  (used  for  decomposition)  was  14  taps.  In  order  to 
increase  the  resolution  in  the  frequency  domain,  we  needed  to  interpolate  between  the 
original  samples.  This  was  done  by  padding  zeros  in  filter  /(n),  and  then  computing  its 
Fourier  transform.  Because  this  filter  was  not  symmetric  with  respect  to  the  origin,  we 
compensated  with  a  |  sample  phase  shift.  The  final  filters  used  in  the  study  are  shown  in 
Table  2.4.3. 


2.4.4  Experimental  Results 

To  date,  our  algorithm  has  been  tested  on  several  examples.  Figures  28(a)  and  28(b)  show 
an  original  signal  and  noisy  signal,  respectively.  The  SNR  is  17.6  dB.  The  denoised  signal 
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h( 

n) 

g( 

k( 

n) 

n 

Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 

-6 

-0.003384 

-0.001523 

0.003384 

-0.001523 

-0.003384 

-0.001523 

-5 

-0.000295 

0.000655 

-0.000295 

-0.000655 

0.000295 

-0.000655 

-4 

0.022730 

0.019115 

-0.022730 

0.019115 

0.022730 

0.019115 

-3 

-0.009154 

0.024154 

-0.009154 

-0.024154 

0.009154 

-0.024154 

-2 

-0.084726 

-0.019306 

0.084726 

-0.019306 

-0.084726 

-0.019306 

-1 

0.094778 

-0.035922 

0.094778 

0.035922 

-0.094778 

0.035922 

0 

0.480049 

0.012827 

-0.480049 

0.012827 

0.480049 

0.012827 

1 

0.480049 

0.012827 

0.480049 

-0.012827 

-0.480049 

-0.012827 

2 

0.094778 

-0.035922 

-0.094778 

-0.035922 

0.094778 

-0.035922 

3 

-0.084726 

-0.019306 

-0.084726 

0.019306 

0.084726 

0.019306 

4 

-0.009154 

0.024154 

0.009154 

0.024154 

-0.009154 

0.024154 

5 

0.022730 

0.019115 

0.022730 

-0.019115 

-0.022730 

-0.019115 

6 

-0.000295 

0.000655 

0.000295 

0.000655 

-0.000295 

0.000655 

7 

-0.003384 

-0.001523 

-0.003384 

0.001523 

0.003384 

0.001523 

l(n) 

n 

Real 

Imaginary 

0 

0.744441 

-0.009214 

1 

0.152932 

0.010146 

2 

0.003917 

0.025936 

3 

-0.038336 

0.001404 

4 

-0.002509 

-0.022426 

5 

0.013544 

-0.014009 

l(-n) 

=  l(n),  n  = 

1,  •••,  13. 

6 

0.002328 

0.001338 

7 

-0.003416 

0.002981 

8 

-0.001008 

-0.000301 

9 

0.000320 

-0.000616 

10 

0.000049 

0.000059 

11 

-0.000048 

0.000099 

12 

0.000002 

0.000002 

13 

0.000005 

-0.000005 

Table  4:  Coefficients  of  h(n),  g(n),  k(n),  and  l(n). 
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is  shown  in  Figure  29(c).  The  new  SNR  is  30.8  dB.  Figures  29(a)  and  29(b)  show  an 
arbitrary  cross  section  of  the  RMI  phantom  image  containing  a  mass.  We  successfully 
preserved  the  mass  while  reducing  noise.  A  two-dimensional  example  is  shown  in  Figure 
30(a)  which  shows  the  original  phantom  image.  It  is  difficult  to  see  the  mass  in  this  image. 
However,  after  denoising,  we  clearly  see  two  masses  on  the  upper  left  and  lower  right 
corners  of  the  image  as  shown  in  Figure  30(b). 

2.4.5  Summary 

Linear  phase  filters  are  very  important  for  mammographic  image  processing.  Such  filters 
enable  control  over  phase  shifts  during  processing.  By  virtue  of  symmetry  (antisymmetry), 
complex  Daubechies  wavelets  were  designed  to  serve  this  advantage  in  analysis. 

In  this  study,  we  proposed  a  denoising  algorithm  based  on  complex  Daubechies 
wavelets.  This  algorithm  aimed  at  removing  noise  while  preserving  features.  We  designed 
an  easy  method  to  classify  waveform  types  and  replaced  the  wavelet  coefficients  of  features 
with  “ideal”  waveform  basis.  The  dc  components  were  also  modified  using  the  unique 
properties  of  complex  Daubechies  wavelets.  Our  approach  differs  from  previous  methods 
including  Mallat’s  [24,  13]  which  used  local  extrema  and  reconstructed  a  denoised  signal 
via  alternative  convex  projection  iteratively.  Despite  high  redundancy  of  frame 
representation,  our  algorithm  was  implemented  efficiently. 

2.5  Segmentation  of  Masses  Using  a  Continuous  Scale 
Representation 

2.5.1  Introduction 

In  this  chapter  of  the  report,  we  describe  a  wavelet  transform  where  a  continuous  scale 
analysis  is  carried  out  through  discretization  and  adaptively  finds  the  scale  at  which  a 
signal  (mass)  is  best  represented. 

Edges  are  an  important  feature  for  analyzing  the  properties  of  images.  Psychophysical 
experiments  suggest  that  the  human  visual  cortex  functions  as  a  population  of  feature 
detectors  tuned  to  edges  and  bars  of  various  widths  and  orientations  [58].  It  has  also  been 
shown  that  visual  information  is  processed  in  parallel  by  a  number  of 
spatial-frequency-tuned  channels  [59]. 

Based  on  such  physiological  findings,  Marr  [60]  developed  an  edge  detection  scheme 
based  on  multiscale  analysis  using  a  Laplacian  of  Gaussian  filter.  Shortly  afterward, 

Witkin  [61]  decomposed  a  signal  into  a  tree  representation  in  terms  of  extrema  of  the  signal 
or  its  derivatives.  This  representation  was  obtained  by  scale-space  filtering  and  it  is 
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complete  under  certain  assumptions. 

In  our  study,  we  have  focused  on  multiscale  wavelet  analysis.  The  wavelet  kernel  is 
predominantly  localized  in  a  certain  region  of  the  space-frequency  plane.  What 
distinguishes  a  wavelet  transform  from  a  short-time  Fourier  transform  is  that  a  wavelet 
may  have  a  larger  spatial  window  for  low  frequency  analysis  and  a  smaller  spatial  window 
for  high  frequency  analysis.  This  property  enables  us  to  zoom  in  on  sharp  variations  and 
design  multiscale  image  processing  algorithms  tuned  for  digital  mammography.  Mallat  first 
developed  a  recursive  algorithm  to  decompose  an  image  into  dyadic  scales  and  reconstruct 
an  original  image  from  its  scale  space  representations  [56].  He  argued  that  when  using  a 
Gaussian  as  a  smoothing  function,  detecting  the  zero-crossings  of  a  wavelet  transform  is 
equivalent  to  a  Marr-Hildreth  edge  detector  [60]. 

A  dyadic  wavelet  analysis  provides  a  complete  non-redundant  representation  [56]. 
However,  here  we  suggest  that  for  the  analysis  of  digital  mammograms,  it  is  desirable  to 
carry  out  analysis  across  arbitrary  scales,  especially,  when  a  certain  mammographic  feature 
lies  between  two  dyadic  scales.  Unser  and  Aldroubi  recently  proposed  a  way  to  calculate  an 
arbitrary  integer  scale  wavelet  transform  using  B-splines  to  approximate  a  wavelet  [62].  In 
this  report,  we  derive  a  scheme  to  calculate  the  wavelet  transform  of  a  discrete  signal  at  an 
arbitrary  scale. 

Although  wavelet  transforms  have  proven  to  be  a  powerful  tool  for  the  enhancement  of 
digital  mammography  [9],  little  is  known  about  choosing  the  optimal  set  of  parameters  for 
a  specific  algorithm  [63].  For  enhancement  of  digital  mammography,  one  must  choose  a 
scale,  threshold  value,  and  gain  factor  in  order  to  target  a  specific  feature  [64] .  To  optimize 
a  wavelet  algorithm,  one  has  to  systematically  investigate  the  effect  of  the  level,  gain,  and 
threshold  [64].  In  addition,  a  distinct  set  of  optimal  parameters  may  be  necessary  for  each 
digital  radiograph. 

Many  image  processing  and  computer  vision  algorithms  are  surprisingly  fragile  when  it 
comes  to  choosing  various  parameters  and  thresholds.  There  have  been  no  good  solutions 
to  setting  the  proper  value  of  the  scale  parameter  in  regularization  networks.  Local 
adaptation  is  needed  in  these  networks  to  make  them  less  brittle.  Using  an  assumption  of 
zero-mean,  additive  Gaussian  noise,  it  can  be  shown  that  the  scale  parameter  for  edge 
detection  problems  should  be  set  to  be  proportional  to  the  variance  of  the  additive  noise 
[65].  However,  these  assumptions  are  too  restrictive  for  most  purposes,  including  digital 
mammography.  More  promising  methods  from  Bayesian  analysis  are  under  investigation  for 
setting  the  scale  and  other  parameters  by  standard  regularization-based  methods  [66,  67]. 

The  Laplacian  filter  is  a  popular  edge  detector  in  image  processing.  However, 
knowledge  of  the  image  feature  size  is  necessary  when  choosing  the  proper  scale  for  feature 
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detection.  In  addition,  traditional  edge  detectors  require  a  smoothing  prefilter  and 
choosing  the  proper  amount  of  smoothing  may  be  difficult.  Compared  to  our  previous 
work,  traditional  image  enhancement  algorithms  may  strengthen  features  within  certain 
spatial  scales  but  choosing  an  optimal  spatial  scale  is  not  obvious. 

It  is  common  to  make  some  assumptions  about  feature  characteristics.  For  example,  we 
can  assume  that  microcalcifications  in  mammograms  occur  within  finer  scales  since  these 
features  are  characterized  by  small,  high-contrast  spots.  Unfortunately,  it  is  more  difficult 
to  make  such  general  assumptions  for  masses  in  mammograms  since  masses  may  assume 
distinct  shapes  and  sizes,  e.g.,  spiculated,  round,  or  irregular  shapes.  Thus,  it  is  necessary 
to  find  a  way  to  identify  an  “optimal”  scale  for  representation  of  masses.  Once  detected, 
their  regions  can  be  marked  and  provide  an  area  of  local  support  for  multiscale  feature 
enhancement,  including  reinforcement  of  the  “halo”  effect,  to  assist  in  the  visibility  of 
masses. 

In  this  chapter,  we  describe  a  method  to  segment  a  mass  from  its  background  using  a 
continuous  multiscale  analysis.  First,  a  suspicious  region  is  identified,  which  may  or  may 
not  contain  a  lesion.  Next,  the  region  is  used  as  a  matched  filter  to  select  a  wavelet  basis. 

A  soft  threshold  is  then  applied.  If  the  selected  region  contains  a  feature,  the 
mammographic  feature  is  segmented  from  its  background  for  further  processing. 


2.5.2  Traditional  Wavelet  Analysis 


A  wavelet  transform  is  a  linear  operator  that  projects  a  signal  onto  a  set  of  basis  functions 
(wavelets).  Each  wavelet  is  a  dilated  and  translated  version  of  a  mother  wavelet.  A  mother 
wavelet  is  a  function  G  such  that 


/OO 

ip{t)dt  =  0. 

-OO 


The  dilation  of  xj^{t)  by  a  factor  a  can  be  described  by 

=  i’l’  (D  ■ 

The  convolution  operator  is  defined  by 

/OO 

f{u)g{t  -  u)du. 

-OO 

The  wavelet  transform  of  a  function  f{t)  at  scale  a  and  position  b  is  defined  as 


Wam=f*Mi)\i=b- 


The  continuous  wavelet  transform  satisfies  energy  conservation  and  f{t)  can  be  perfectly 
reconstructed  from  its  wavelet  transform  coefficients  [68].  Both  scale  a  and  position  b  may 
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vary  continuously  over  the  set  of  real  numbers.  At  a  fine  scale  uq?  fbe  support  of  the  basis 
V’ao(^)  is  small  and  wavelet  coefficients  Waof{b)  are  sensitive  to  detail.  For  a  larger  scale  Ci, 
the  support  of  V’ai(0  has  greater  extent  and  wavelet  coefficients  Wa^fib)  are  more  sensitive 
to  slow  changing  information  (or  larger  objects). 

For  a  particular  class  of  wavelets,  the  scale  parameter  can  be  sampled  on  a  dyadic  grid 
without  losing  the  mathematical  completeness  property  of  the  transform.  A  dyadic  wavelet 
transform  of  a  continuous  function  f(t)  at  scale  2^  and  at  position  b  is  defined  as  [13]. 

W2kf{b)  =  f*M-t)\t=b.  (85) 


A  function  f{t)  can  then  be  reconstructed  from  its  transform  coefficients  alone  by: 

OO  OO 

f{t)=  E  W2uf*X2'>{t)=  E  (86) 

k-^oo  k=-oo 

where  wavelets  and  x(0  satisfy  the  condition 

CO 

^  ^(2'=a;)A’(2M  =  1,  (8^) 

k=—oo 

and  ^(w)  and  X{(^)  are  the  Fourier  transforms  of  and  x(0?  respectively. 

The  dyadic  wavelet  transform  filters  used  in  this  report  were  designed  analytically  as 
follows. 

Similar  to  [13],  we  define  a  smoothing  function  $(w)  such  that  [25] 

+  00 

$(a;,)  =  n  (88) 

p  =  l 


Since  Ut=i  cos(2-Pa>)  =  we  let  [10,  25] 


$(cc;) 

H{u;) 


and 


(89) 

(90) 


Choosing  (?{co)  =  (4isin(|))^  =  -16sin^(|),  the  wavelet  il>(u’)  can  be  defined  simply  as  [25] 


<Hu)  =  G 


(91) 


Filter  K{^)  is  then  defined  by 


K(u) 


l-|-g(u-)P 

a(w) 


2^ 

16 


(l  +  cos 


(92) 
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Finally,  the  reconstructing  wavelet  y(t)  can  be  formulated  as 


=  (93) 

We  note  that  since  0(t<;)  =  jg  a  cubic  spline  whose  integral  is  equal  to  1,  is  a 

4 

second  order  derivative  of  a  cubic  spline  smoothing  function.  For  more  detailed  treatment 
of  this  class  of  wavelets  please  refer  to  Section  2.2  of  this  report. 

2.5.3  An  Arbitrary  Scale  Algorithm  for  the  Discrete  Case 

A  variation  of  Mallat’s  algorithm  for  discrete  dyadic  analysis  is  presented  next.  Mallat’s 
algorithm  can  analyze  a  discrete  signal  x[n],  of  finite  energy  [13].  In  his  algorithm,  the 
scaling  function  $(a;)  must  satisfy  the  constraint 

+  00 

yuj  e  R,  0  <  Cl  <  ^  |$(ci;  +  2n7r)p  <  where  Ci  and  C2  are  constants.  (94) 

n=z—oo 

Marr  [60]  claimed  also  that  there  exists  some  function  f{t)  €  L’^{R)  (not  unique)  such  that 
[13] 

'■  ■'  roo 

\/n  e  Z,  Sifin)  =  /  -  n)dt  =  x[n\.  (95) 

J—oo 

At  each  scale  2-^,  we  may  decompose  S^jf[k]  into  S^j+ifik]  =  bj^k  and  W^j+iflk]  =  Cj^k  [13]. 
The  decomposition  algorithm  is: 
j  =  0 

while  {j  <  J) 

=  Si  f*G, 

SU^f  =  Sif*H, 
i  =  i  + 1 

end  of  while, 

where  Hp{uj)  =  H{2^u))  and  Gp{u!)  =  G{2^u;).  The  reconstruction  algorithm  is: 

j  =  J 

while  (j  >  0) 

si-j  =  wif  *  Kj.i  +  sif  *  i?,-i 

i  =  i  - 1 

end  of  while, 

where  Kp{oj)  =  K{2'^lo). 
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In  frequency  domain,  this  can  be  formulated  as: 


(96) 

(97) 


=  F{iv)^u) 

=  S^JG{2^u) 

=  Sl,fH{2^cv) 

=  Si-^f  H{2^-'^u;)H{2^u^) 

=  F{io)^io  )  n 

2=0 

=  F(u;)^{2oj  )  n 

2  =  1 

=  F{u)^{2^+^uj)  (98) 


where  the  frequency  response  of  the  analysis  filter  0{u)  is 


(101) 


Notice  that  in  the  algorithm  above,  filters  at  scale  2^  are  obtained  by  putting  2^  ^  —  1 
zeros  between  consecutive  coefficients  of  filters  at  level  1.  However,  this  introduces 
undesirable  high  frequency  components  into  the  analysis  filter  0[lij)  as  shown  in  Figure 
31(a).  Thus,  a  modification  was  introduced  as  shown  in  Figure  31(b).  By  including  an 
ideal  lowpass  filter  with  cutoff  frequency  7r/2-^-\  a  mathematically  perfect  reconstruction  is 
sacrified  to  obtain  superior  frequency  performance  as  shown  in  Figure  31(c).  We  can  show 
that  perfect  reconstruction  is  almost  achievable  with  respect  to  contrast  resolution 
currently  provided  by  digital  mammography.  When  using  the  wavelet  defined  previously, 
the  corresponding  analysis  filter  for  our  algorithm  becomes 


^7n(2^'+'  LO) 


'F(2-^+^  uj)  when  \u}  +  2kT^\  <  2^ 

0  when  \u  +  2k'K\  >  2^~^'k. 


(102) 
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Frequency  bands  for  differerrt  levels 


Frequency  bands  for  different  levels 


(a)  (c) 


(b) 


Figure  31;  (a)  The  frequency  bands  for  standard  analysis,  (b)  Block  diagram  of  our  modified 
algorithm  (3  levels  shown),  (c)  The  frequency  response  for  the  subbands  of  our  algorithm. 


If  satisfies  the  admissibility  condition,  $(w)  =  0  when  a;  =  0.  Thus  —  0  when 

a;  =  0.  And  so  we  can  define 


0(2^0;)  = 


$(w)  ' 


(103) 


This  analysis  filter  and  its  frequency  response  are  comparable  to  our  original  analysis 
filters  as  shown  in  Figure  32.  Our  filter  appears  more  smooth,  which  is  desirable  for  image 


10  20  30  40 

(d)  spatial 


Figure  32:  (a)  Standard  analysis  filter’s  frequency  response  at  level  3.  (b)  Our  analysis 
filter’s  frequency  response  at  level  3.  (c)  Second  derivative  wavelet  at  level  3.  (d)  Our  new 
wavelet  at  level  3. 


analysis.  In  addition,  the  ripples  caused  by  the  ideal  low  pass  filter  are  very  subtle.  Our 
new  algorithm  also  enables  us  to  zoom  in-between  octave  levels  and  trace  a  particular 
feature  in  a  mammogram.  Let  us  now  define 


^(2aoj)  when  |a;  +  2A:7r|  <  27r/a, 
0  when  |a;  -|-  2A:7r|  >  27r/a, 


(104) 


and 


D{aLo) 


^m(2aa;) 

$(w)  ’ 


where  a  >  2  or  a  =  1,  when  a  =  2-’,  D{au)  =  0{2^lv)  as  shown  in  Figure  33. 


(105) 


2.5,4  Locating  the  Best  Scale 

In  the  previous  section  we  showed  how  wavelet  coefficients  can  be  calculated  at  an  arbitrary 
scale.  We  next  discuss  a  method  for  determining  an  appropriate  scale  based  on  local 
feature  sizes  within  an  image.  Using  the  image  as  a  matched  filter,  finite  number  of  scales 
are  searched  to  find  the  wavelet  coefficients  that  most  resemble  the  features  of  interest. 
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dOGhdot  line:  a*6.6667 


Figure  33:  Capability  to  zoom  in  between  2  levels. 
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Figure  34:  (a)  The  best  scale  to  detect  a  bump  in  a  simulated  mass  profile,  (b)  The  best 
scale  to  detect  a  bump  in  a  noisy  image  profile. 


This  method  is  demostrated  in  Figure  34.  The  left-hand  side  of  Figure  34(a)  shows 
three  scales  of  an  image  profile,  containing  bumps  of  distinct  width.  We  calculated  the 
wavelet  transform  for  each  image.  At  each  scale,  the  maximum  of  wavelet  coefficients 
across  the  shifting  parameter  was  found.  The  relation  between  scale  and  the  wavelet 
maximum  is  shown  on  the  right-hand  side  of  Figure  34(a).  Note  that  there  is  exactly  one 
scale  a*  corresponding  to  the  maximum  of  wavelet  coefficients  across  both  shifting  and 
scale  parameters.  The  value  a*  is  clearly  a  function  of  the  size  of  the  feature  within  the 
cropped  region. 

Figure  34(b)  shows  an  image  profile  corrupted  with  random  noise.  Our  method  is  still 
able  to  identify  candidates  for  the  best  scale  in  the  presence  of  severe  noise. 

2.5.5  Application  of  the  Segmentation  Algorithm 

X-ray  images  of  an  RMI  model  156  phantom  (Radiation  Measurements  Inc.,  Middleton, 
WI)  were  investigated.  The  RMI  phantom  contains  five  masses  and  six  fibrils  which  mimic 
two  objects  of  interest  in  mammography.  The  radiographic  image  quality  of  this  phantom, 
as  well  as  the  total  number  of  objects  which  are  deemed  to  be  visible,  form  a  key 
component  of  the  mammography  accreditation  program  administered  by  the  American 
College  of  Radiology  (ACR).  The  RMI  phantom  is  designed  so  that  at  least  one  object  in 
each  category  is  not  visible  and  one  object  is  at  the  borderline  of  the  visibility  threshold 
[69]. 

An  X-ray  image  of  the  RMI  phantom  is  ideal  for  studying  image  segmentation  since  it 
has  well  defined  features  of  clinical  interest  which  are  present  both  above  and  below  the 
threshold  for  visibility.  Figure  35(a)  shows  a  schematic  representation  of  the  insert  in  the 
RMI.  For  imaging  purposes,  this  insert  is  placed  in  block  of  acrylic  designed  to  simulate  an 
average  compressed  breast.  Figure  35(b)  shows  an  example  of  the  corresponding 
radiographic  image  obtained  with  a  conventional  mammography  unit  as  used  for  ACR 
accreditation. 

We  suggest  that  traditional  multiscale  analysis  carried  out  at  dyadic  scales  is  flawed  in 
the  sense  of  best  representation  and  sensitivity.  To  demonstrate  graphically  the  flaw  of 
processing  at  dyadic  scales  alone  we  present  the  following  seqence  of  images,  obtained  by 
processing  a  digital  radiograph  of  the  RMI  phantom  (shown  in  Figure  35)  containing  three 
masses  (small,  medium,  and  large).  Figure  36(a)  shows  regions  obtained  by  the 
zero-crossings  of  wavelet  coefficients  extracted  from  the  dyadic  scale  five.  Note  that  the 
bounded  regions  corresponding  to  the  large  and  medium  masses  are  quite  deformed  and 
there  are  many  small  ’’false  positives”  detected  at  this  level  of  analysis.  However,  at  the 
next  dyadic  scale  (level  six)  shown  in  Figure  36(b),  the  boundaries  of  the  small  and  large 
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Figure  35:  RMI  phantom:  (a)  internal  schematic;  (b)  digital  radiograph. 
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Visibility  on  x-ray  image 

Eeature 

Mass  thickness/  diameter 

Fibril  (nylon  fiber  diameter) 

Easy  (E) 

750  fim  /  6  mm 

750  nm 

Borderline  (B) 

500  nm  /  5  mm 

540  yum 

Invisible  (I) 

250  nm  /  3  mm 

400  yum 

Table  5:  Characteristics  of  masses  and  fibrils  investigated. 


masses  are  merged  with  an  adjacent  L-shaped  border  detected  within  the  phantom. 

Thus,  a  better  representation  may  live  (exist)  in  between  the  two  dyadic  scales,  level 
five  (which  proved  to  be  too  fine  an  analysis)  and  six  (too  coarse  an  analysis).  Indeed,  this 
is  exactly  what  we  find.  Figure  37(a)  shows  the  regions  corresponding  to  zero-crossings 
extracted  at  a  scale  inbetween  levels  five  and  six.  Note  the  shape  of  each  of  the  three 
masses  are  distinct  and  well  defined.  Figure  37(b),  shows  the  three  masses  detected  by 
applying  simple  geometric  constraints  to  each  labelled  segmented  region  within  the  best 
level  of  analysis.  The  one  false  positive  present  is  most  likely  due  to  the  lack  of 
sophisitication  of  our  present  (simple)  algorithm. 

Figure  35(a)  shows  that  the  phantom  also  contains  objects  which  simulate 
microcalcifications,  however  this  study  focused  on  fibrils  and  masses.  Table  5  summarizes 
the  key  parameters  of  the  masses  and  fibrils  that  were  investigated.  Identification  of  each 
object  is  by  category  of  fibril  (F)  or  mass  (M)  and  by  its  visibility  of  easy  (E),  borderline 
(B),  or  invisible  (I). 

Radiographic  images  of  the  RMI  phantom  were  obtained  using  a  Digital  Spot 
Mammography  (DSM)  system  attached  to  a  Lorad  Breast  Biopsy  System  (Lorad 
Corporation,  Danbury,  CT).  The  detector  area  was  50  mm  by  50  mm  which  is  smaller  than 
the  RMI  phantom  size  shown  in  Figure  35.  As  a  result,  only  20  %  was  able  to  be  captured 
in  any  single  digital  image.  The  DSM  image  matrix  size  used  in  this  study  was  512  by  512 
so  that  the  nominal  pixel  size  was  100  /^m.  Separate  images  were  obtained  of  mass  and 
fibril  regions  of  the  phantom.  For  the  mass  region,  technique  factors  of  22  kVp  and  16  mAs 
were  used,  whereas  for  the  fibril  region,  22  kVp  and  80  mAs  were  used.  Additional  images 
were  also  acquired  at  32  kVp  and  200  mAs,  where  the  increased  radiation  reduced  the 
noise  and  permitted  the  B/I  masses  and  fibrils  to  become  visible.  By  this  process,  the 
locations  of  the  invisible  and  borderline  masses  in  test  images  (standard  radiation  dose) 
were  established  with  a  high  level  of  accuracy. 

The  acquired  image  data  set  was  transferred  to  a  SUN  workstation  (Mountain  View, 
Ca.)  for  subsequent  analysis  and  display.  A  256  by  256  pixel  region  of  interest 
incorporating  the  E  and  B  masses  was  extracted  for  subsequent  processing  whereas  for  the 
invisible  mass,  a  cropped  region  of  183  by  183  was  used.  For  the  fibrils,  the  selected  region 
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Figure  36:  (a)  Blobs  detected  from  dyadic  wavelet  coefficients  at  level  5,  (b)  Blobs  detected 
from  dyadic  wavelet  coefficients  at  level  6. 


Figure  37:  (a)  Suspicious  regions  detected  for  wavelet  coefficients  at  scale  28.5.  (b)  Suspi¬ 
cious  region  identified  as  possible  masses  providing  local  support  for  enhancement  of  wavelet 
coefficients. 


was  128  by  128  pixels.  The  choice  of  these  regions  was  determined  by  the  need  to  eliminate 
all  extraneous  objects  for  this  part  of  the  image  processing  exercise.  For  comparison 
purposes,  a  uniform  background  of  200  by  200  pixels  containing  no  objects  was  also 
investigated,  and  categorized  as  none  (N)  below. 

Segmentation  -  Mass 

Figure  38  shows  representative  images  of  the  cropped  regions  corresponding  to  three 
masses  (E,  B  and  I)  as  well  as  a  background  region  with  no  objects  (N).  These  images  were 
generated  using  22  kVp  and  16  mAs.  Also  shown  are  the  same  four  images  after  histogram 
equalization  and  unsharp  masking  enhancement.  As  expected,  the  processed  images  in 
Figure  35  shows  negligible  improvement  in  overall  mass  visibility  using  either  of  these  two 
common  image  enhancement  techniques. 

Figure  39  shows  the  segmentation  results  obtained  for  an  easily  visible  mass,  M-E.  A 
profile  through  the  center  of  M-E  is  depicted  in  Eigure  39(a),  and  the  corresponding  profile 
of  the  wavelet  coefficients  is  shown  in  Figure  40(b).  It  is  the  interscale  value,  a  =  48,  which 
best  defines  the  wavelet  coefficient  profile  and  therefore  best  segments  the  mass  (see  Figure 
39(d)). 

Figure  41  shows  the  segmentation  results  obtained  for  a  mass  with  borderline  visibility, 
M  —  B.  A  profile  through  the  center  of  M  —  B  is  depicted  in  Figure  42(a),  and  the 
corresponding  profile  of  the  wavelet  coefficient  is  shown  in  Figure  42(b).  Note  that  the 
interscale  value  a  =  28  best  defines  the  wavelet  coefficient  profile  and  therefore  best 
segments  the  mass  (see  Figure  41(e)).  The  segmentation  algorithm  clearly  defined  this 
mass. 

Figure  43  shows  the  segmentation  results  obtain  for  the  invisible  mass,  M  —  I.  A  profile 
through  the  center  of  M  -  /  is  depicted  in  Figure  44(a),  and  the  corresponding  profile  of 
the  wavelet  coefficients  is  shown  in  Figure  44(b).  Here  the  interscale  value  a  =  28  again 
best  defines  the  wavelet  coefficient  profile  and  therefore  best  segments  the  mass  (see  Figure 
43(e)).  It  is  clear  that  the  segmentation  algorithm  developed  in  this  work  was  capable  of 
finding  the  mass  which  is  normally  invisible  on  standard  radiographs  of  this  phantom. 

Figure  45  shows  the  results  obtained  when  processing  the  background  region  with  no 
object  in  the  field  of  view. 

Segmentation  -  Fibrils 

Figure  47  shows  representative  images  of  the  cropped  regions  corresponding  to  three 
fibrils  (E,  B,  and  I).  As  in  the  case  of  masses,  these  images  were  generated  using  34  kVp 
(22  is  preferred)  and  56  mAs  (16  mAs  is  preferred).  Also  shown  are  these  same  images 
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Figure  38:  Different  processing  methods,  window  and  leveling  (W/L),  histogram  equalization 
(HE),  unsharp  masking  (UM).  (a)  The  easy  case,  (b)  The  borderline  case,  (c)  The  invisible 


Profile  of  a  mass  (large) 


Figure  40:  Profiles  of  the  large  mass  and  its  wavelet  coefficients 


Profile  of  a  mass  (small) 


Figure  44:  Profiles  of  the  small  mass  and  its  wavelet  coefficients 


Profile  of  a  mass  (small) 


Figure  46:  Profiles  of  the  small  mass  and  its  wavelet  coefficients. 


using  histogram  equalization  as  well  as  an  unsharp  masking  enhancement.  The  processed 
images  in  Figure  47  show  that  little  improvement  in  overall  fibril  visibility  is  achieved  using 
either  of  these  three  standard  image  enhancement  techniques. 

Figure  48  shows  the  segmentation  results  obtained  for  an  easily  visible  fibril,  F-E.  A 
profile  through  the  center  of  F-E  is  depicted  in  Figure  49(a),  and  the  corresponding  profile 
of  wavelet  coefficients  is  shown  in  Figure  49(b).  It  is  the  interscale  value  a  =  7.5  which  best 
defines  the  wavelet  coefficient  profile  and  therefore  best  segments  the  shape  of  the  fibril 
(see  Figure  48(a)). 

Figure  50  shows  the  segmentation  results  obtained  for  fibril  with  borderline  visibility, 
F-B.  A  profile  through  the  center  of  F-B  is  depicted  in  Figure  51(a),  and  the  corresponding 
profile  of  the  wavelet  coefficients  is  shown  in  Figure  51(b).  The  interscale  value  a  =  10  best 
defined  the  wavelet  coefficient  profile  and  therefore  best  captured  the  fibril  (see  Figure 
50(e)).  The  segmentation  algorithm  clearly  defined  this  fibril. 

Figure  52  shows  the  segmentation  results  obtained  for  a  fibril  that  is  invisible,  F-I.  A 
profile  through  the  center  of  F-I  is  depicted  in  Figure  53(a),  and  the  corresponding  profile 
of  the  wavelet  coefficients  is  shown  in  Figure  53(b).  It  is  an  interscale  value  which  best 
defines  the  wavelet  coefficient  profile  and  therefore  best  captures  the  fibril  (see  Figure 
52(a)).  The  segmentation  algorithm  clearly  defined  this  subtle  fibril. 

2.5.6  Summary 

We  have  shown  that  regions  corresponding  to  masses  can  be  identified  reliably  through 
frame  representations  of  a  continuous  multiscale  analysis.  We  showed  that  subtle  features 
characteristic  of  mammographic  findings  required  a  finer  parameterization  of  scale  space 
than  provided  by  traditional  methods  of  wavelet  analysis. 

Basis  functions  for  carrying  out  continuous  multiscale  analysis  were  designed  to  be 
symmetric  and  have  zero-phase  providing  closed  contours  (via  zero- crossings)  of  emergent 
features  within  each  level  of  scale. 

Spatial  locality  was  preserved  via  a  frame  representation  enabling  segmentation  of 
masses  and  adaptive  contrast  enhancement.  We  suggest  that  masses  located  in  dense  and 
somewhat  dense  mammograms  can  be  reliably  characterized  and  identified  by  applying 
geometric  constraints  (size  and  shape)  to  the  set  of  segmented  regions. 
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Profile  of  a  fibril  (easy  a) 


20 


Figure  50:  Fibril  (B-borderline)  and  its  edges:  (a)  original  image;  (b)-(d)  dyadic  scales;  (e) 
edges  at  best  scale. 


Profile  of  a  fibril  (hard  b) 


Figure  53:  Profiles  of  the  fibril  and  its  wavelet  coefficients. 


2.6  Enhancement  of  Computer  Simulated  Masses  for 
Mammography  Using  Wavelet  Analysis 

2.6.1  Introduction 

In  this  section  we  report  results  on  the  application  of  enhancement  by  hard  thresholding. 
We  present: 

•  mathematical  details  of  modeling  of  image  formation, 

•  details  of  hard  thresholding  method,  and 

•  methods  (mathematical  and  psychophysical)  for  evaluation  of  the  enhancement 
technique. 

In  section  2.6.3  we  present  the  enhancement  evaluation  methods.  In  section  2.6.4  we 
demonstrate  results  of  our  experimentation  with  the  hard  thresholding  technique  and  in 
section  2.6.5  we  present  conclusions  and  future  directions. 

2.6.2  Methodology 

Phantoms  Simulated  by  Mathematical  Construction 

The  phantom  used  in  this  study  was  designed  to  replicate  circularly  shaped  masses  in  a 
uniform  noisy  background. 

1.  Image:  For  computational  purposes,  a  512  by  512  matrix  containing  a  number  of 
simulated  masses  was  used.  There  were  8  bits  per  pixel  corresponding  to  256  grey 
levels.  The  phantom  was  divided  into  25  squares  (a  5  by  5  grid).  The  size  of  each 
square  was  100  by  100  pixels.  Thirteen  of  these  squares  were  chosen  randomly  to 
contain  a  simulated  mass  while  the  remaining  twelve  contained  only  noise. 

2.  The  masses:  Each  mass  was  modelled  as  a  circle  with  a  radius  equal  to  r.  The  center 
of  the  mass  coincided  with  the  center  of  the  square.  The  signal  intensity  of  the  mass 
was  computed  by: 

I{x,  y)  =  a%J (r2  -  -  y'^)  (106) 

where  x  and  y  were  the  coordinates  with  respect  to  a  local  coordinate  system  (point 
of  origin  for  the  local  coordinate  system  was  the  geometrical  center  of  the  mass).  The 
parameter  a  controls  the  local  contrast  and  r  controls  the  size  of  the  simulated  mass 
(r  is  the  radius  of  the  mass).  The  background  had  a  grey  level  of  value  128. 


99 


PSF  vs  normalized  distance 


Figure  54:  A  plot  of  the  Point  Spread  Function  (PSF)  as  a  function  of  the  normalized 
distance  d/r^. 

3.  Noise  and  blurring:  In  mammography,  the  dominant  source  of  image  noise  is 

quantum  mottle  and  for  this  simulation,  a  lOmi?  exposure  level  at  the  detector  was 
used  with  an  average  photon  energy  of  20  keV .  A  42  nm  pixel  size  was  used  since 
this  corresponds  to  a  Ak  x  6k  matrix  size  for  a  conventional  25  cm  x  30  cm 
mammogram  [70].  As  a  result  the  mean  number  of  photons  per  pixel  should  be  a 
1000  which  corresponds  to  a  standard  deviation  in  the  noise  level  of  3%  [71]).  The 
limiting  resolution  in  screen-film  mammography  is  in  the  range  15  to  20  line  pairs  per 
mm  [72].  The  computer  simulated  phantom  was  convolved  with  a  smoothing  function 
to  simulate  a  reduced  modulation  transfer  function.  A  screen  thickness  of  80  f^m  was 
assumed  and  a  two  dimensional  Point  Spread  Function  (PSF)  generated  using  an 
expression  derived  by  Barrett  and  Swindell  [73].  A  plot  of  the  Point  Spread  Function 
used  in  our  simulations  as  a  function  of  the  normalized  distance  rjds-,  is  presented  in 
Figure  54. 

Figure  55  (a)  shows  one  of  twenty-five  squares  that  was  cropped  from  a  simulated 
phantom.  The  size  of  the  square  is  100  by  100  pixels.  This  square  contains,  located  at  its 
center,  a  mass  of  radius  of  20  pixels  and  it  is  corrupted  with  zero  mean  Gaussian  noise  of 
standard  deviation  of  15.  Figures  55  (b)  and  55  (c)  show  two  phantoms  of  the  same  size 
with  only  a  mass  and  only  noise.  Figures  56  (a),  56  (b)  and  56  (c)  show  three  profiles 
chosen  from  the  Figure  55  (a).  Finally,  Figure  57  shows  a  typical  simulated  phantom 
containing  thirteen  masses. 
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Figure  57:  A  typical  mathematical  phantom  used  in  our  study.  The  size  of  the  phantom 
was  512  by  512  pixels.  The  phantom  contained  13  masses  and  was  corrupted  with  random 
Gaussian  noise. 
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Function  usgcJ  for  the  thresholding  of  the  wavelet  coefficients 


Figure  58:  A  plot  of  the  piecewise  linear  function  of  the  wavelet  output  /{w)  as  a  function 
of  the  input  w  of  the  wavelet  coefficients. 

Image  Enhancement 

Thresholding  of  the  wavelet  coefficients  for  the  purpose  of  denoising  is  a  popular  theme 
in  current  literature.  In  this  work  we  investigated  the  use  of  hard  thresholding  for  the 
denoising  of  the  simulated  phantoms.  In  this  method  for  each  decomposition  level  i  a 
threshold  value  U  was  selected  and  all  the  coefficients  below  that  level  were  set  to  zero. 

The  rest  remained  unaltered.  Mathematically  this  can  be  expressed  as  follows: 

^  f  0  for  <  wi,i  <  tl 
^  \  'Wi,i  otherwise, 

where  wi^i  is  the  /-th  wavelet  coefficient  of  the  i-th  level  t”  and  are  the  positive  and 
negative  thresholds  respectively  for  decomposition  level  i.  Figure  58  shows  schematically 
the  thresholding  strategy.  The  choice  of  the  thresholds  and  was  based  on  the 
percentage  of  the  wavelet  coefficients  that  remain  unaltered  after  the  thresholding.  In  other 
words  if  we  set  the  threshold  at  20  %,  and  were  chosen  in  such  a  way  as  to  leave  80  % 
of  the  negative  wavelet  coefficients  and  80  %  percent  of  the  positive  wavelet  coefficients 
unaltered.  The  equality: 

f.  =  -t-  (108) 

does  not  necessarily  hold  since  the  histogram  of  the  coefficients  for  all  levels  is  not 
symmetric  with  respect  to  the  ?/  =  0  axis.  Figures  59(a)-(f)  show  the  histograms  of  wavelet 
coefficients  of  the  six  decomposition  levels  for  the  phantom  of  Figure  57.  The  vertical  lines 
denote  the  thresholds  for  the  percentages  0  %,  10  ,  100  %. 

After  performing  the  inverse  transform  the  resulting  image  was  window  leveled  for 
better  visualization  and  contrast. 
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Figure  60:  An  image  I  is  defined  in  the  rectangular  domain  D.  A  mass  m  is  defined  in  the 
domain  Dm  (white  region).  The  area  surrounding  the  domain  Dm,  D  —  Dm,  is  used  for  the 
calculation  of  the  background  B  (black  region) . 

2.6.3  Evaluation  of  Enhancement 

In  this  Section  we  describe  our  methodology  for  quantitative  and  qualitative  evaluation.  A 
discussion  of  the  relation  between  the  Enhancement  Factor  EF,  and  the  performance  of  an 
observer  has  been  presented  in  a  previous  publication  [63]. 


Enhancement  Factor 

The  mean  signal  contrast,  C ,  in  each  simulated  mass  was  defined  as: 

^  _  IPm  _  Ip-Dm  (109) 

IPm  Id -Dm 

where  S  is  the  original  signal  and  dA  is  the  area  differential.  Dm  is  the  area  where  the 
signal  (mass)  resides,  D  —  Dm  is  the  area  that  surrounds  the  mass  and  D  is  the  area  of  the 
square.  For  the  calculation  of  C  we  used  the  following  formula  [74]: 


_ ^ ^  (110) 

13  ^  \7rr^  100^  —  Trr^/ 

where  Sj  is  the  gross  count  in  the  j-th  square  in  the  area  where  the  signal  resided  {Dm) 
and  Bj  was  the  gross  count  in  the  same  square  where  the  Sj  was  calculated  but  for  the 
area  that  surrounds  the  simulated  mass  (see  Figure  60).  Equation  110  is  a  discretized 
version  of  Equation  109.  The  summation  considers  the  blocks  that  contain  a  mass. 

An  estimation  of  the  noise  N,  in  the  image  is  provided  by  the  following  formula: 


7V  = 


\ 


12 


12 

E 

i=i 


(111) 


106 


where  aj  is  the  standard  deviation  in  the  square  where  no  signal  resides.  The  summation 
for  that  case  considers  the  blocks  that  contain  only  noise. 

These  definitions  for  contrast  and  noise  were  used  to  obtain  the  input  contrast  to  noise 
ratio  (CNRi)  and  the  corresponding  value  of  the  output  contrast  to  noise  ratio  {CNRo)  in 
the  processed  image.  The  resultant  enhancement  factor,  EF,  is  given  by  the  expression: 

CNR,  _  {C/N\ 

CNRi  {ClN)i  ■ 


Psychophysics  Evaluations 

The  performance  of  the  mass  enhancement  algorithm  was  evaluated  through 
psychophysical  experiments.  The  phychophysical  technique  used  for  the  evaluation  of  the 
improvement  of  the  visibility  was  the  two- alternative  forced  choice  (2 AFC)  measurement. 
Recent  results  suggest  that  the  alternative  forced  choice  (AFC)  techniques  are  superior 
when  the  experimental  conditions  can  be  precisely  controled  using  synthetic 
computer- generated  images  or  images  produced  with  phantoms  and  diagnostic  equipment 
[75].  The  2 AFC  method  was  adopted,  instead  of  the  4AFC  (which  according  to  Burgess 
[75]  is  the  best  choice),  because  of  its  simplicity  and  its  reduced  computational  and  storage 
costs. 

In  this  method  a  radiologist  was  given  two  images  of  size  128  by  128  pixels.  The  noise 
in  both  images  was  normal  and  of  the  same  variance.  One  of  the  two  images  was  assigned  a 
mass  of  known  intensity,  shape  and  size.  The  two  images  had  equal  probabilities  of  being 
assigned  the  mass.  A  reference  image,  indicating  the  correct  position,  shape  and  strength 
of  the  mass  in  a  noiseless  background  was  present  to  assist  the  radiologist  in  his  task  of 
locating  the  mass. 

The  same  process  was  repeated  for  both  the  original  and  processed  images.  Eight  point 
CNRs  were  selected  and  for  each  point  100  images  were  reviewed.  The  radiologist  used  a 
mouse  to  indicate  his  preference  on  the  right  or  left  image.  For  each  point  the  proportion 
of  the  correct  responses,  was  calculated  and  plotted  as  a  function  of  the  CNR.  For  the 
evaluation  of  the  processed  image  we  generated  images,  one  with  a  mass  and  noise  and  the 
other  with  only  noise.  Both  images  were  processed  using  the  threshold  algorithm  and 
presented  as  the  two  alternatives  to  the  radiologist. 

Figure  61  shows  the  graphical  user  interface  built  for  the  2AFC  experiments.  The  top 
row  contains  the  two  images  (left  and  right)  the  radiologist  is  forced  to  choose  from.  In  the 
middle  row  (on  the  right  side)  one  can  see  the  reference  image  containing  the  noiseless 
mass  (reference  image).  This  image  refreshes  the  radiologist’s  memory  about  the  shape, 
size  and  position  of  the  simulated  mass.  On  the  left  corner  there  is  the  control  panel.  It 
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contains  the  buttons  the  radiologist  has  to  click  on  in  order  to  choose  the  left  or  the  right 
image,  a  counter  showing  the  number  of  the  image  and  a  button  he/ she  can  use  if  he/ she 
wants  to  review  a  previous  image. 

2.6.4  Results 

Choosing  Optimized  Thresholds 

Our  task  was  to  choose  the  threshold  values  for  each  level  in  such  a  way  that  we 
obtained  best  possible  improvement  in  visibility.  As  a  measure  of  improvement  we  used  the 
E F  described  earlier.  In  order  to  obtain  this  optimum  choice  we  relied  on  an  empirical 
study  of  the  variation  of  the  EF  with  the  choice  of  the  threshold  value  for  each  level. 

In  order  to  reduce  the  computational  complexity  we  grouped  the  six  decomposition 
levels  into  three  groups: 

Group  1  (high  decomposition  levels):  Levels  1  and  2.  This  group  contained  noise  and  the 
fine  details  of  the  image. 

Group  2  (medium  decomposition  levels):  Levels  3  and  4.  This  group  contained  a  mixture 
of  noise  and  signal. 

Group  3  (low  decomposition  levels):  Levels  5  and  6.  This  group  contained  mainly  pure 
signal. 

Our  study  was  performed  in  two  stages.  First  we  studied  the  effects  of  the  threshold 
levels  of  each  group  and  then  we  concentrated  on  the  areas  and  levels  we  found  most 
important  and  performed  a  more  detailed  study.  In  the  first  part  of  our  study  eleven 
threshold  levels  were  chosen,  thresholding  at  0  %,  10  %,. . . ,  100  %  percent  of  the 
coefficients  for  that  level.  The  EF  was  plotted  as  a  mesh  as  a  function  of  the  thresholds  of 
two  groups.  This  was  performed  for  all  combinations.  Results  from  the  first  part  of  our 
study  are  presented  in  Figure  62.  These  results  suggested  that  only  a  very  small  percentage 
of  the  wavelet  coefficients  of  the  high  levels  should  not  be  thresholded.  Therefore  we 
concentrated  in  the  area  of  threshold  values  from  90  %  to  100  %.  Figure  63  shows  the  plot 
of  the  EF  as  a  function  of  the  threshold  values  for  Groups  1  and  3.  It  is  obvious  that  even 
if  1  %  of  the  wavelet  coefficients  survives  the  thresholding  there  is  no  significant 
improvement  on  the  EF.  We  therefore  concluded  that  all  the  wavelet  coefficients  of  Levels  1 
and  2  should  be  set  to  zero. 

For  groups  2  and  3  we  pursued  into  a  more  detailed  study  to  obtain  an  estimation  of 
the  optimized  threshold  values.  For  group  2  we  concentrated  on  the  area  from  90  %  to 
100  %  and  found  that  an  optimum  choice  depends  heavily  on  the  CNR.  We  therefore  opted 
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Figure  61:  The  graphical  user  interface  used  during  the  2 AFC  experiments. 


EF  as  a  function  of  the  thresholds  of  Groups  1  and  3 


Figure  63:  Plot  of  the  EF  as  a  function  of  the  thresholds  of  Groups  1  and  3.  For  Group  1 
we  focused  on  the  area  from  91  %  to  99  %.  Notice  that  the  increase  of  the  EF  is  minimal. 
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for  a  conservative  choice  of  the  threshold  and  set  99  %  of  the  coefficients  of  levels  3  and  4 
equal  to  zero.  After  obtaining  the  optimized  threshold  values  for  Groups  1  and  2  we 
studied  the  effects  of  the  choice  of  the  threshold  values  of  levels  5  and  6  seperately.  Figure 
64  (a)  shows  a  plot  of  the  EF  as  a  function  of  the  threshold  values  of  Levels  5  and  6. 
Figures  64  (b)— (c)  shows  profiles  of  the  same  plot.  In  summary  we  found  that  the  optimum 
choice  of  threshold  values  were:  level  1:  100  %,  level  2:  100  %,  level  3:  99  %,  level  4:  99  %, 
level  5:  87  %,  level  6:  75  %.  We  point  out  that  as  it  is  obvious  from  Figure  64  that  the 
improvement  in  CNR  does  not  show  a  big  variation  with  the  choice  of  threshold  value  for 
Levels  5  and  6.  Even  if  the  threshold  value  for  that  level  is  not  close  to  optimum,  one 
would  obtain  a  significant  improvement  in  CNR. 

Results  on  Simulated  Masses 

In  this  subsection  the  improvement  of  the  visibility  of  masses  in  synthetic  phantoms 
using  the  thresholding  algorithm  are  presented.  In  Figure  65a  an  128  by  128  pixel  image 
containing  a  mass  is  shown.  Figures  65  (c)-(h)  show  the  coefficients  that  resulted  from  the 
wavelet  decomposition  of  that  image.  The  image  presented  in  that  figure  was  cropped  from 
a  512  by  512  pixel  simulated  phantom.  The  image  of  the  wavelet  coefficients  have  been 
window- leveled  for  best  display.  The  radius  of  the  mass  was  equal  to  15  pixels  and  the 
noise  had  a  standard  deviation  equal  to  10.  The  processed  image  is  shown  in  Figure  65  (b). 
The  resulting  CNR  was  found  to  be  equal  to  28  (EF=39).  Another  example  is  presented  in 
Figure  66  (a).  The  CNR  of  this  image  was  0.1,  and  was  also  cropped  from  a  512  by  512 
pixel  image.  Figure  66  (b)  shows  the  enhanced  image.  It  is  obvious  that  even  in  this 
extreme  case  there  is  a  significant  improvement  in  the  visibility  of  the  mass.  For  this  case 
the  EF  was  found  to  be  equal  to  25. 

Psychophysics  Experiments 

The  2AFC  method  was  used  to  evaluate  the  enhancement  of  the  visibility  of  masses. 
The  resulting  True  Positive  Curves  for  a  radiologist  are  plotted  as  function  of  the  CNR  (see 
Figure  67).  We  note  that  both  True  Positive  Curves  have  a  sigmoid  shape  starting  from 
100  %,  for  low  CNRs  (i.e.,  the  radiologist  is  able  to  see  everything)  and  ending,  for  high 
CNRs,  to  50  %  (i.e.,  the  radiologist  has  no  clue  as  to  where  the  mass  is  and  he  randomly 
picks  an  image).  We  can  clearly  see  that  the  transition  region  for  all  three  radiologists  has 
been  shifted  to  higher  CNRs.  This  indicates  a  real  improvement  in  the  visibility  of  the 
masses. 
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EF  vs.  thresholds  of  5th  and  6th  level 


EF  as  a  function  of  the  threshold  of  the  5th  level 
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EF  as  a  function  of  the  threshold  of  the  6th  level 
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Figure  64;  (a)  A  plot  of  the  EF  as  a  function  of  the  thresholds  of  Level  5  and  Level  6;  (b)  a 
plot  of  the  EF  as  a  function  of  the  threshold  of  level  5  (profile  of  (a));  (c)  a  plot  of  the  EF 
as  a  function  of  the  threshold  of  level  6  (profile  of  (a)). 


Figure  65:  (a)  Original  image,  (b)  Processed  image,  (c)-(h)  Six  levels  of  wavelet  decompo¬ 
sition. 
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(b) 

Figure  66:  (a)  A  simulated  phantom  containing  a  mass.  The  CNR  was  0.1.  (b)  The  resulting 
image  after  the  enhancement.  The  CNR  was  equal  to  25. 
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2.6.5  Summary 

Substantial  progress  was  made  in  the  development  of  techniques  to  assist  in  the 
quantitative  evaluation  of  wavelet  based  image  processing  algorithms.  We  showed  that 
synthetic  images  can  provide  a  useful  tool  towards  benchmarking  the  effectiveness  of 
multiscale  enhancement  algorithms.  However  the  conclusions  reached  about  the 
effectiveness  of  these  methods  largely  depend  on  the  realism  of  the  synthetic  data.  Our 
method  capitalized  on  our  previous  results  to  introduce  a  more  realistic  model  of  image 
formation.  During  the  next  year  we  plan  to  incorporate  a  model  for  the  appearance  of 
tissue,  such  as  presented  by  Karssemeijer  in  his  work  for  calcification  detection  in  digital 
mammograms. 

Finally,  we  reported  on  the  development  of  objective  ways  to  assess  the  performance  of 
wavelet  image  processing  algorithms.  An  important  consideration  in  our  enhancement 
algorithms  was  the  selection  of  threshold  values.  We  studied  methods  for  estimation  of  an 
optimum  threshold  value  and  carried  out  an  empirical  threshold  selection  scheme.  Our 
experiments  showed  that  by  deleting  wavelet  coefficients  within  the  high  frequency 
channels  (e.g.,  levels  1  and  2)  we  achieved  considerable  improvement  of  the  contrast  to 
noise  ratio.  Aggressively  thresholding  wavelet  coefficients  within  the  lower  spatial 
frequency  channels  had  the  effect  of  removing  some  of  the  features  of  the  image  (e.g.,  edges 
which  appear  exactly  within  these  levels).  It  is  therefore  concluded  that  a  more  clever 
scheme  is  needed  that  will  not  threshold  wavelet  coefficients  with  respect  to  spatial  locality 
but  will  use  spatial/frequency  information  when  modifying  coefficients  for  denoising  and 
enhancement. 
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3  Conclusions 


During  the  past  year,  we  have  made  significant  progress  in  the  development  of  a 
methodology  for  accomplishing  adaptive  contrast  enhancement  by  multiscale 
representations.  Our  studies  have  demonstrated  that  features  extracted  from 
multiresolution  representations  can  provide  an  adaptive  mechanism  for  the  local  emphasis 
of  salient  and  subtle  features  of  importance  to  mammography.  The  improved  contrast  of 
mammographic  features  make  these  techniques  appealing  for  computed  aided  diagnosis 
(CAD)  and  screening  mammography.  Screening  mammography  examinations  are  certain  to 
grow  substantially  in  the  next  few  years,  and  analytic  methods  that  can  assist  general 
radiologists  in  reading  mammograms  shall  be  of  great  importance. 

In  the  paragraphs  below,  we  summarize  our  progress  and  identify  future  directions  of 
research  to  be  carried  out  during  the  next  year  of  our  investigation. 

We  have  described  some  underlying  criteria  to  be  followed  when  designing  a  database  of 
digitized  mammograms.  Respecting  these  criteria  we  constructed  a  database  of  105 
abnormal  cases  (420  images)  with  a  total  of  118  biopsy  proven  masses.  All  cases  were 
supplied  by  Shands  Hospital  at  the  University  of  Florida,  an  ACR  accredited  and  FDA 
certified  institution.  Mammograms  were  digitized  with  a  Lumiscan  75  scanner  (Sunnyvale, 
California)  at  116-micron  pixel  resolution,  12-bit  gray  scale  resolution  and  optical  density 
range  of  0.05-3.00.  All  abnormalities  were  annotated  hierarchically  by  a  radiologist  from 
Shands  Hospital  experienced  in  reading  mammograms  using  an  inhouse  tool  called  XMaiti 
described  previously  in  our  1994  annual  report.  Annotations  include  the  location,  shape, 
margins,  density,  size,  surrounding  glandular  pattern,  mammography  finding  and  diagnosis 
based  on  the  standardized  American  College  of  Radiology  Breast  Imaging  Reporting  and 
Data  System.  In  addition,  we  developed  a  subjective  evaluation  of  the  detectability  of  a 
lesion  and  evaluated  all  lesions  based  on  this  scale.  The  current  database  consists  of  78 
malignant  masses  (66.1  %)  and  40  benign  masses  (33.9  %)  distributed  according  to 
subtlety  as  follows:  28.8  %  obvious,  28  %  relatively  obvious,  25.4  %  subtle,  14.4  %  very 
subtle,  and  3.4  %  extremely  subtle. 

We  computed  the  coefficients  for  filters  of  low  order  splines  used  as  smoothing  functions 
in  a  fast  hierarchical  digital  filtering  implementation  of  a  discrete  dyadic  wavelet  transform 
for  both  one  and  two-dimensional  processing.  A  mirror  extended  input  signal  to  a  filter 
bank  implementation  of  the  discrete  dyadic  wavelet  transform  enabled  us  to  take  advantage 
of  the  symmetry/antisymmetry  of  both  signals  and  filters  for  an  efficient  implementation  of 
a  dyadic  transform.  There  are  no  restrictions  on  the  size  of  the  input  signal  which  makes 
this  scheme  attractive  for  extracting  arbitrary  rectangular  regions  from  digital 
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mammograms.  The  presented  algorithm  is  iterative  across  scales:  a  noniterative  scheme 
suitable  for  high  speed  (interactive)  parallel  processing  can  be  derived  easily  by  using  the 
described  approach  on  each  dyadic  scale  independently. 

We  have  shown  the  limitation  of  traditional  discrete  wavelet  transforms  (DWT)  for 
characterizing  band-limited  features,  including  subtle  findings  in  mammographic  images. 

To  more  reliably  isolate  noise  and  identify  features  through  scale  space,  we  formulated  and 
developed  sub-octave  wavelet  transforms.  If  the  number  of  sub-octave  bands  in  each  octave 
is  a  power  of  2,  we  showed  that  the  sub-octave  wavelet  transform  of  a  function  can  be 
implemented  by  a  single  pair  of  analysis  and  synthesis  wavelets.  FIR  filters  for  a  class  of 
wavelets  designed  for  a  DWT  can  be  used  to  implement  sub-octave  wavelet  analysis.  A 
method  of  wavelet  shrinkage  for  noise  suppression  was  presented.  An  algorithm  for  the 
generalized  case  of  adaptive  gain  processing  was  developed  and  shown  to  enhance  features 
without  amplifying  noise.  These  efforts  have  significance  in  supporting  interactive 
processing  while  yielding  increased  sensitivity  of  mammographic  findings. 

Our  future  efforts  will  include  developing  feature  clustering  techniques  and  carrying  out 
statistical  analysis  of  multiscale  coefficients  for  more  efficient  removal  of  noise  and 
enhancement  of  features.  We  will  also  explore  iterative  methods  for  identifying  coherent 
structures  of  masses  (selected  from  our  mammographic  database)  that  are  persistent 
through  scale  space  representations  derived  from  our  method  of  sub-octave  wavelet  analysis. 

We  have  always  advocated  the  use  of  linear  phase  filters  for  mammographic  image 
processing  and  analysis.  Such  filters  enable  control  over  phase  shifts  and  limit  the 
introduction  of  artifacts  during  non-linear  processing.  By  virtue  of  their  symmetry 
(antisymmetry),  complex  Daubechies  wavelets  were  designed  to  serve  this  advantage  in 
image  analysis. 

We  described  a  denoising  algorithm  based  on  complex  Daubechies  wavelets.  The 
algorithm  was  aimed  at  removing  noise  while  preserving  mammographic  features.  We 
designed  a  simple  method  to  classify  waveform  types  and  replaced  wavelet  coefficients  of 
features  with  (mathematically)  ideal  waveform  basis.  The  DC  components  of  the  analysis 
were  modified  using  some  unique  properties  of  complex  Daubechies  wavelets. 

We  have  shown  that  regions  corresponding  to  masses  can  be  identified  reliably  through 
frame  representations  of  a  continuous  multiscale  analysis.  We  showed  that  subtle  features 
characteristic  of  mammographic  findings  required  a  finer  parameterization  of  scale  space 
than  provided  by  traditional  methods  of  wavelet  analysis  carried  out  at  dyadic  scales. 

Basis  functions  for  carrying  out  continuous  multiscale  analysis  were  designed  to  be 
symmetric  and  have  zero-phase  providing  closed  contour  boundaries  (via  zero-crossings)  of 
emergent  features  within  each  level  of  scale. 
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spatial  locality  was  preserved  via  a  frame  representation  enabling  segmentation  of 
masses  and  adaptive  contrast  enhancement.  We  suggest  that  masses  in  dense 
mammograms  can  be  reliably  characterized  and  identified  via  geometric  constraints  (size 
and  shape)  within  a  best  scale.  This  shall  be  pursued  during  the  next  year  of  research. 

Substantial  progress  was  made  in  the  development  of  techniques  to  assist  in  the 
quantitative  evaluation  of  wavelet  based  image  processing  algorithms.  We  showed  that 
synthetic  images  can  provide  a  useful  tool  towards  benchmarking  the  effectiveness  of 
multiscale  enhancement  algorithms.  However  the  conclusions  reached  about  the 
effectiveness  of  these  methods  largely  depend  on  the  realism  of  the  synthetic  data.  Our 
method  capitalized  on  our  previous  results  to  introduce  a  more  realistic  model  of  image 
formation.  During  the  next  year  we  plan  to  incorporate  a  model  for  the  appearance  of 
tissue,  such  as  presented  by  Karssemeijer  in  his  work  for  calcification  detection  in  digital 
mammograms. 

Finally,  we  reported  on  the  development  of  objective  ways  to  assess  the  performance  of 
wavelet  image  processing  algorithms.  An  important  consideration  in  our  enhancement 
algorithms  was  the  selection  of  threshold  values.  We  studied  methods  for  estimation  of  an 
optimum  threshold  value  and  carried  out  an  empirical  threshold  selection  scheme.  Our 
experiments  showed  that  by  deleting  wavelet  coefficients  within  the  high  frequency  channels 
(e.g.,  levels  1  and  2)  we  achieved  considerable  improvement  of  the  contrast  to  noise  ratio. 
Aggressively  thresholding  wavelet  coefficients  within  the  lower  spatial  frequency  channels 
had  the  effect  of  removing  some  of  the  features  of  the  image  (e.g.,  edges  which  appear 
exactly  within  these  levels).  It  is  therefore  concluded  that  a  more  clever  scheme  is  needed 
that  will  not  threshold  wavelet  coefficients  with  respect  to  spatial  locality  but  will  use 
spatial/frequency  information  when  modifying  coefficients  for  denoising  and  enhancement. 

These  results  are  encouraging  and  continue  to  support  that  wavelet  based  image 
processing  algorithms  can  play  an  important  role  in  improving  the  imaging  performance  of 
digital  mammography.  In  summary,  we  have  exceeded  the  goals  as  described  in  our  Task 
III  and  Task  IV  Statement  of  Work  for  and  our  research  and  development  plans  remain  on 
schedule. 
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